ONEWAY: Fix printf format specifier.
[pspp] / src / language / stats / kruskal-wallis.c
1 /* Pspp - a program for statistical analysis.
2    Copyright (C) 2010, 2011 Free Software Foundation, Inc.
3
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.
8
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.
13
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/>. */
16
17
18 #include <config.h>
19
20 #include "kruskal-wallis.h"
21
22 #include <gsl/gsl_cdf.h>
23 #include <math.h>
24
25 #include "data/casereader.h"
26 #include "data/casewriter.h"
27 #include "data/dataset.h"
28 #include "data/dictionary.h"
29 #include "data/format.h"
30 #include "data/subcase.h"
31 #include "data/variable.h"
32 #include "libpspp/assertion.h"
33 #include "libpspp/hmap.h"
34 #include "libpspp/message.h"
35 #include "libpspp/misc.h"
36 #include "math/sort.h"
37 #include "output/tab.h"
38
39 #include "gl/minmax.h"
40 #include "gl/xalloc.h"
41
42
43 /* Returns true iff the independent variable lies in the range [nst->val1, nst->val2] */
44 static bool
45 include_func (const struct ccase *c, void *aux)
46 {
47   const struct n_sample_test *nst = aux;
48
49   if (0 < value_compare_3way (&nst->val1, case_data (c, nst->indep_var), var_get_width (nst->indep_var)))
50     return false;
51
52   if (0 > value_compare_3way (&nst->val2, case_data (c, nst->indep_var), var_get_width (nst->indep_var)))
53     return false;
54
55   return true;
56 }
57
58
59 struct rank_entry
60 {
61   struct hmap_node node;
62   union value group;
63
64   double sum_of_ranks;
65   double n;
66 };
67
68 /* Return the entry with the key GROUP or null if there is no such entry */
69 static struct rank_entry *
70 find_rank_entry (const struct hmap *map, const union value *group, size_t width)
71 {
72   struct rank_entry *re = NULL;
73   size_t hash  = value_hash (group, width, 0);
74
75   HMAP_FOR_EACH_WITH_HASH (re, struct rank_entry, node, hash, map)
76     {
77       if (0 == value_compare_3way (group, &re->group, width))
78         return re;
79     }
80   
81   return re;
82 }
83
84 /* Calculates the adjustment necessary for tie compensation */
85 static void
86 distinct_callback (double v UNUSED, casenumber t, double w UNUSED, void *aux)
87 {
88   double *tiebreaker = aux;
89
90   *tiebreaker += pow3 (t) - t;
91 }
92
93
94 struct kw
95 {
96   struct hmap map;
97   double h;
98 };
99
100 static void show_ranks_box (const struct n_sample_test *nst, const struct kw *kw, int n_groups);
101 static void show_sig_box (const struct n_sample_test *nst, const struct kw *kw);
102
103 void
104 kruskal_wallis_execute (const struct dataset *ds,
105                         struct casereader *input,
106                         enum mv_class exclude,
107                         const struct npar_test *test,
108                         bool exact UNUSED,
109                         double timer UNUSED)
110 {
111   int i;
112   struct ccase *c;
113   bool warn = true;
114   const struct dictionary *dict = dataset_dict (ds);
115   const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test, parent);
116   const struct caseproto *proto ;
117   size_t rank_idx ;
118
119   int total_n_groups = 0.0;
120
121   struct kw *kw = xcalloc (nst->n_vars, sizeof *kw);
122
123   /* If the independent variable is missing, then we ignore the case */
124   input = casereader_create_filter_missing (input, 
125                                             &nst->indep_var, 1,
126                                             exclude,
127                                             NULL, NULL);
128
129   input = casereader_create_filter_weight (input, dict, &warn, NULL);
130
131   /* Remove all those cases which are outside the range (val1, val2) */
132   input = casereader_create_filter_func (input, include_func, NULL, 
133         CONST_CAST (struct n_sample_test *, nst), NULL);
134
135   proto = casereader_get_proto (input);
136   rank_idx = caseproto_get_n_widths (proto);
137
138   /* Rank cases by the v value */
139   for (i = 0; i < nst->n_vars; ++i)
140     {
141       double tiebreaker = 0.0;
142       bool warn = true;
143       enum rank_error rerr = 0;
144       struct casereader *rr;
145       struct casereader *r = casereader_clone (input);
146
147       r = sort_execute_1var (r, nst->vars[i]);
148
149       /* Ignore missings in the test variable */
150       r = casereader_create_filter_missing (r, &nst->vars[i], 1,
151                                             exclude,
152                                             NULL, NULL);
153
154       rr = casereader_create_append_rank (r, 
155                                           nst->vars[i],
156                                           dict_get_weight (dict),
157                                           &rerr,
158                                           distinct_callback, &tiebreaker);
159
160       hmap_init (&kw[i].map);
161       for (; (c = casereader_read (rr)); case_unref (c))
162         {
163           const union value *group = case_data (c, nst->indep_var);
164           const size_t group_var_width = var_get_width (nst->indep_var);
165           struct rank_entry *rank = find_rank_entry (&kw[i].map, group, group_var_width); 
166
167           if ( NULL == rank)
168             {
169               rank = xzalloc (sizeof *rank);
170               value_clone (&rank->group, group, group_var_width);
171
172               hmap_insert (&kw[i].map, &rank->node,
173                            value_hash (&rank->group, group_var_width, 0));
174             }
175
176           rank->sum_of_ranks += case_data_idx (c, rank_idx)->f;
177           rank->n += dict_get_case_weight (dict, c, &warn);
178
179           /* If this assertion fires, then either the data wasn't sorted or some other
180              problem occured */
181           assert (rerr == 0);
182         }
183
184       casereader_destroy (rr);
185
186       /* Calculate the value of h */
187       {
188         struct rank_entry *mre;
189         double n = 0.0;
190
191         HMAP_FOR_EACH (mre, struct rank_entry, node, &kw[i].map)
192           {
193             kw[i].h += pow2 (mre->sum_of_ranks) / mre->n;
194             n += mre->n;
195
196             total_n_groups ++;
197           }
198         kw[i].h *= 12 / (n * ( n + 1));
199         kw[i].h -= 3 * (n + 1) ; 
200
201         kw[i].h /= 1 - tiebreaker/ (pow3 (n) - n);
202       }
203     }
204
205   casereader_destroy (input);
206   
207   show_ranks_box (nst, kw, total_n_groups);
208   show_sig_box (nst, kw);
209
210   /* Cleanup allocated memory */
211   for (i = 0 ; i < nst->n_vars; ++i)
212     {
213       struct rank_entry *mre, *next;
214       HMAP_FOR_EACH_SAFE (mre, next, struct rank_entry, node, &kw[i].map)
215         {
216           hmap_delete (&kw[i].map, &mre->node);
217           free (mre);
218         }
219       hmap_destroy (&kw[i].map);
220     }
221
222   free (kw);
223 }
224
225 \f
226 #include "gettext.h"
227 #define _(msgid) gettext (msgid)
228
229
230 static void
231 show_ranks_box (const struct n_sample_test *nst, const struct kw *kw, int n_groups)
232 {
233   int row;
234   int i;
235   const int row_headers = 2;
236   const int column_headers = 1;
237   struct tab_table *table =
238     tab_create (row_headers + 2, column_headers + n_groups + nst->n_vars);
239
240   tab_headers (table, row_headers, 0, column_headers, 0);
241
242   tab_title (table, _("Ranks"));
243
244   /* Vertical lines inside the box */
245   tab_box (table, 1, 0, -1, TAL_1,
246            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
247
248   /* Box around the table */
249   tab_box (table, TAL_2, TAL_2, -1, -1,
250            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
251
252   tab_text (table, 1, 0, TAT_TITLE, 
253             var_to_string (nst->indep_var)
254             );
255
256   tab_text (table, 3, 0, 0, _("Mean Rank"));
257   tab_text (table, 2, 0, 0, _("N"));
258
259   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
260   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
261
262
263   row = column_headers;
264   for (i = 0 ; i < nst->n_vars ; ++i)
265     {
266       int tot = 0;
267       const struct rank_entry *re;
268
269       if (i > 0)
270         tab_hline (table, TAL_1, 0, tab_nc (table) -1, row);
271       
272       tab_text (table,  0, row,
273                 TAT_TITLE, var_to_string (nst->vars[i]));
274
275       HMAP_FOR_EACH (re, const struct rank_entry, node, &kw[i].map)
276         {
277           struct string str;
278           ds_init_empty (&str);
279
280           var_append_value_name (nst->indep_var, &re->group, &str);
281
282           tab_text   (table, 1, row, TAB_LEFT, ds_cstr (&str));
283           tab_double (table, 2, row, TAB_LEFT, re->n, &F_8_0);
284           tab_double (table, 3, row, TAB_LEFT, re->sum_of_ranks / re->n, 0);
285
286           tot += re->n;
287           row++;
288           ds_destroy (&str);
289         }
290       tab_double (table, 2, row, TAB_LEFT,
291                   tot, &F_8_0);
292       tab_text (table, 1, row++, TAB_LEFT, _("Total"));
293     }
294
295   tab_submit (table);
296 }
297
298
299 static void
300 show_sig_box (const struct n_sample_test *nst, const struct kw *kw)
301 {
302   int i;
303   const int row_headers = 1;
304   const int column_headers = 1;
305   struct tab_table *table =
306     tab_create (row_headers + nst->n_vars * 2, column_headers + 3);
307
308   tab_headers (table, row_headers, 0, column_headers, 0);
309
310   tab_title (table, _("Test Statistics"));
311
312   tab_text (table,  0, column_headers,
313             TAT_TITLE | TAB_LEFT , _("Chi-Square"));
314
315   tab_text (table,  0, 1 + column_headers,
316             TAT_TITLE | TAB_LEFT, _("df"));
317
318   tab_text (table,  0, 2 + column_headers,
319             TAT_TITLE | TAB_LEFT, _("Asymp. Sig."));
320
321   /* Box around the table */
322   tab_box (table, TAL_2, TAL_2, -1, -1,
323            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
324
325
326   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
327   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
328
329   for (i = 0 ; i < nst->n_vars; ++i)
330     {
331       const double df = hmap_count (&kw[i].map) - 1;
332       tab_text (table, column_headers + 1 + i, 0, TAT_TITLE, 
333                 var_to_string (nst->vars[i])
334                 );
335
336       tab_double (table, column_headers + 1 + i, 1, 0,
337                   kw[i].h, 0);
338
339       tab_double (table, column_headers + 1 + i, 2, 0,
340                   df, &F_8_0);
341
342       tab_double (table, column_headers + 1 + i, 3, 0,
343                   gsl_cdf_chisq_Q (kw[i].h, df),
344                   0);
345     }
346
347   tab_submit (table);
348 }