1 /* Pspp - a program for statistical analysis.
2 Copyright (C) 2010 Free Software Foundation, Inc.
4 This program is free software: you can redistribute it and/or modify
5 it under the terms of the GNU General Public License as published by
6 the Free Software Foundation, either version 3 of the License, or
7 (at your option) any later version.
9 This program is distributed in the hope that it will be useful,
10 but WITHOUT ANY WARRANTY; without even the implied warranty of
11 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12 GNU General Public License for more details.
14 You should have received a copy of the GNU General Public License
15 along with this program. If not, see <http://www.gnu.org/licenses/>. */
20 #include "kruskal-wallis.h"
22 #include <gsl/gsl_cdf.h>
25 #include <data/casereader.h>
26 #include <data/casewriter.h>
27 #include <data/dictionary.h>
28 #include <data/format.h>
29 #include <data/procedure.h>
30 #include <data/subcase.h>
31 #include <data/variable.h>
33 #include <libpspp/assertion.h>
34 #include <libpspp/message.h>
35 #include <libpspp/misc.h>
36 #include <libpspp/hmap.h>
37 #include <math/sort.h>
45 include_func (const struct ccase *c, void *aux)
47 const struct n_sample_test *nst = aux;
49 if (0 < value_compare_3way (&nst->val1, case_data (c, nst->indep_var), var_get_width (nst->indep_var)))
52 if (0 > value_compare_3way (&nst->val2, case_data (c, nst->indep_var), var_get_width (nst->indep_var)))
61 struct hmap_node node;
68 static struct rank_entry *
69 find_rank_entry (const struct hmap *map, const union value *group, size_t width)
71 struct rank_entry *re = NULL;
72 size_t hash = value_hash (group, width, 0);
74 HMAP_FOR_EACH_WITH_HASH (re, struct rank_entry, node, hash, map)
76 if (0 == value_compare_3way (group, &re->group, width))
84 distinct_callback (double v UNUSED, casenumber t, double w UNUSED, void *aux)
86 double *tiebreaker = aux;
88 *tiebreaker += pow3 (t) - t;
98 static void show_ranks_box (const struct n_sample_test *nst, const struct kw *kw, int n_groups);
99 static void show_sig_box (const struct n_sample_test *nst, const struct kw *kw);
102 kruskal_wallis_execute (const struct dataset *ds,
103 struct casereader *input,
104 enum mv_class exclude,
105 const struct npar_test *test,
112 const struct dictionary *dict = dataset_dict (ds);
113 const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test, parent);
114 const struct caseproto *proto ;
117 int total_n_groups = 0.0;
119 struct kw *kw = xcalloc (nst->n_vars, sizeof *kw);
121 /* If the independent variable is missing, then we ignore the case */
122 input = casereader_create_filter_missing (input,
127 input = casereader_create_filter_weight (input, dict, &warn, NULL);
129 /* Remove all those cases which are outside the range (val1, val2) */
130 input = casereader_create_filter_func (input, include_func, NULL,
131 CONST_CAST (struct n_sample_test *, nst), NULL);
133 proto = casereader_get_proto (input);
134 rank_idx = caseproto_get_n_widths (proto);
136 /* Rank cases by the v value */
137 for (i = 0; i < nst->n_vars; ++i)
139 double tiebreaker = 0.0;
141 enum rank_error rerr = 0;
142 struct casereader *rr;
143 struct casereader *r = casereader_clone (input);
145 r = sort_execute_1var (r, nst->vars[i]);
147 /* Ignore missings in the test variable */
148 r = casereader_create_filter_missing (r, &nst->vars[i], 1,
152 rr = casereader_create_append_rank (r,
154 dict_get_weight (dict),
156 distinct_callback, &tiebreaker);
158 hmap_init (&kw[i].map);
159 for (; (c = casereader_read (rr)); case_unref (c))
161 const union value *group = case_data (c, nst->indep_var);
162 const size_t group_var_width = var_get_width (nst->indep_var);
163 struct rank_entry *rank = find_rank_entry (&kw[i].map, group, group_var_width);
167 rank = xzalloc (sizeof *rank);
168 value_clone (&rank->group, group, group_var_width);
170 hmap_insert (&kw[i].map, &rank->node,
171 value_hash (&rank->group, group_var_width, 0));
174 rank->sum_of_ranks += case_data_idx (c, rank_idx)->f;
175 rank->n += dict_get_case_weight (dict, c, &warn);
177 /* If this assertion fires, then either the data wasn't sorted or some other
182 casereader_destroy (rr);
185 struct rank_entry *mre;
188 HMAP_FOR_EACH (mre, struct rank_entry, node, &kw[i].map)
190 kw[i].h += pow2 (mre->sum_of_ranks) / mre->n;
195 kw[i].h *= 12 / (n * ( n + 1));
196 kw[i].h -= 3 * (n + 1) ;
198 kw[i].h /= 1 - tiebreaker/ (pow3 (n) - n);
202 casereader_destroy (input);
204 show_ranks_box (nst, kw, total_n_groups);
205 show_sig_box (nst, kw);
207 for (i = 0 ; i < nst->n_vars; ++i)
209 struct rank_entry *mre, *next;
210 HMAP_FOR_EACH_SAFE (mre, next, struct rank_entry, node, &kw[i].map)
212 hmap_delete (&kw[i].map, &mre->node);
215 hmap_destroy (&kw[i].map);
222 #include <output/tab.h>
224 #define _(msgid) gettext (msgid)
228 show_ranks_box (const struct n_sample_test *nst, const struct kw *kw, int n_groups)
232 const int row_headers = 2;
233 const int column_headers = 1;
234 struct tab_table *table =
235 tab_create (row_headers + 2, column_headers + n_groups + nst->n_vars);
237 tab_headers (table, row_headers, 0, column_headers, 0);
239 tab_title (table, _("Ranks"));
241 /* Vertical lines inside the box */
242 tab_box (table, 1, 0, -1, TAL_1,
243 row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
245 /* Box around the table */
246 tab_box (table, TAL_2, TAL_2, -1, -1,
247 0, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
249 tab_text (table, 1, 0, TAT_TITLE,
250 var_to_string (nst->indep_var)
253 tab_text (table, 3, 0, 0, _("Mean Rank"));
254 tab_text (table, 2, 0, 0, _("N"));
256 tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
257 tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
260 row = column_headers;
261 for (i = 0 ; i < nst->n_vars ; ++i)
264 const struct rank_entry *re;
267 tab_hline (table, TAL_1, 0, tab_nc (table) -1, row);
269 tab_text (table, 0, row,
270 TAT_TITLE, var_to_string (nst->vars[i]));
272 HMAP_FOR_EACH (re, const struct rank_entry, node, &kw[i].map)
275 ds_init_empty (&str);
277 var_append_value_name (nst->indep_var, &re->group, &str);
279 tab_text (table, 1, row, TAB_LEFT, ds_cstr (&str));
280 tab_double (table, 2, row, TAB_LEFT, re->n, &F_8_0);
281 tab_double (table, 3, row, TAB_LEFT, re->sum_of_ranks / re->n, 0);
287 tab_double (table, 2, row, TAB_LEFT,
289 tab_text (table, 1, row++, TAB_LEFT, _("Total"));
297 show_sig_box (const struct n_sample_test *nst, const struct kw *kw)
300 const int row_headers = 1;
301 const int column_headers = 1;
302 struct tab_table *table =
303 tab_create (row_headers + nst->n_vars * 2, column_headers + 3);
305 tab_headers (table, row_headers, 0, column_headers, 0);
307 tab_title (table, _("Test Statistics"));
309 tab_text (table, 0, column_headers,
310 TAT_TITLE | TAB_LEFT , _("Chi-Square"));
312 tab_text (table, 0, 1 + column_headers,
313 TAT_TITLE | TAB_LEFT, _("df"));
315 tab_text (table, 0, 2 + column_headers,
316 TAT_TITLE | TAB_LEFT, _("Asymp. Sig."));
318 /* Box around the table */
319 tab_box (table, TAL_2, TAL_2, -1, -1,
320 0, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
323 tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
324 tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
326 for (i = 0 ; i < nst->n_vars; ++i)
328 const double df = hmap_count (&kw[i].map) - 1;
329 tab_text (table, column_headers + 1 + i, 0, TAT_TITLE,
330 var_to_string (nst->vars[i])
333 tab_double (table, column_headers + 1 + i, 1, 0,
336 tab_double (table, column_headers + 1 + i, 2, 0,
339 tab_double (table, column_headers + 1 + i, 3, 0,
340 gsl_cdf_chisq_Q (kw[i].h, df),