Merge branch 'master' into output
[pspp] / src / language / stats / wilcoxon.c
1 /* Pspp - a program for statistical analysis.
2    Copyright (C) 2008, 2009 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
19 #include <config.h>
20
21 #include "wilcoxon.h"
22
23 #include <gsl/gsl_cdf.h>
24 #include <math.h>
25 #include <signal.h>
26 #include <unistd.h>
27
28 #include <data/casereader.h>
29 #include <data/casewriter.h>
30 #include <data/dictionary.h>
31 #include <data/format.h>
32 #include <data/procedure.h>
33 #include <data/subcase.h>
34 #include <data/variable.h>
35 #include <gsl/gsl_cdf.h>
36 #include <libpspp/assertion.h>
37 #include <libpspp/message.h>
38 #include <libpspp/misc.h>
39 #include <math/sort.h>
40 #include <math/wilcoxon-sig.h>
41 #include <output/table.h>
42 #include <signal.h>
43 #include <unistd.h>
44
45 #include "minmax.h"
46 #include "xalloc.h"
47
48 static double
49 append_difference (const struct ccase *c, casenumber n UNUSED, void *aux)
50 {
51   const variable_pair *vp = aux;
52
53   return case_data (c, (*vp)[0])->f - case_data (c, (*vp)[1])->f;
54 }
55
56 static void show_ranks_box (const struct wilcoxon_state *,
57                             const struct two_sample_test *,
58                             const struct dictionary *);
59
60 static void show_tests_box (const struct wilcoxon_state *,
61                             const struct two_sample_test *,
62                             bool exact, double timer);
63
64
65
66 static void
67 distinct_callback (double v UNUSED, casenumber n, double w UNUSED, void *aux)
68 {
69   struct wilcoxon_state *ws = aux;
70
71   ws->tiebreaker += pow3 (n) - n;
72 }
73
74 #define WEIGHT_IDX 2
75
76 void
77 wilcoxon_execute (const struct dataset *ds,
78                   struct casereader *input,
79                   enum mv_class exclude,
80                   const struct npar_test *test,
81                   bool exact,
82                   double timer)
83 {
84   int i;
85   bool warn = true;
86   const struct dictionary *dict = dataset_dict (ds);
87   const struct two_sample_test *t2s = (struct two_sample_test *) test;
88
89   struct wilcoxon_state *ws = xcalloc (sizeof (*ws), t2s->n_pairs);
90   const struct variable *weight = dict_get_weight (dict);
91   struct variable *weightx = var_create_internal (WEIGHT_IDX, 0);
92   struct caseproto *proto;
93
94   input =
95     casereader_create_filter_weight (input, dict, &warn, NULL);
96
97   proto = caseproto_create ();
98   proto = caseproto_add_width (proto, 0);
99   proto = caseproto_add_width (proto, 0);
100   if (weight != NULL)
101     proto = caseproto_add_width (proto, 0);
102
103   for (i = 0 ; i < t2s->n_pairs; ++i )
104     {
105       struct casereader *r = casereader_clone (input);
106       struct casewriter *writer;
107       struct ccase *c;
108       struct subcase ordering;
109       variable_pair *vp = &t2s->pairs[i];
110
111       ws[i].sign = var_create_internal (0, 0);
112       ws[i].absdiff = var_create_internal (1, 0);
113
114       r = casereader_create_filter_missing (r, *vp, 2,
115                                             exclude,
116                                             NULL, NULL);
117
118       subcase_init_var (&ordering, ws[i].absdiff, SC_ASCEND);
119       writer = sort_create_writer (&ordering, proto);
120       subcase_destroy (&ordering);
121
122       for (; (c = casereader_read (r)) != NULL; case_unref (c))
123         {
124           struct ccase *output = case_create (proto);
125           double d = append_difference (c, 0, vp);
126
127           if (d > 0)
128             {
129               case_data_rw (output, ws[i].sign)->f = 1.0;
130
131             }
132           else if (d < 0)
133             {
134               case_data_rw (output, ws[i].sign)->f = -1.0;
135             }
136           else
137             {
138               double w = 1.0;
139               if (weight)
140                 w = case_data (c, weight)->f;
141
142               /* Central point values should be dropped */
143               ws[i].n_zeros += w;
144               case_unref (output);
145               continue;
146             }
147
148           case_data_rw (output, ws[i].absdiff)->f = fabs (d);
149
150           if (weight)
151            case_data_rw (output, weightx)->f = case_data (c, weight)->f;
152
153           casewriter_write (writer, output);
154         }
155       casereader_destroy (r);
156       ws[i].reader = casewriter_make_reader (writer);
157     }
158   caseproto_unref (proto);
159
160   for (i = 0 ; i < t2s->n_pairs; ++i )
161     {
162       struct casereader *rr ;
163       struct ccase *c;
164       enum rank_error err = 0;
165
166       rr = casereader_create_append_rank (ws[i].reader, ws[i].absdiff,
167                                           weight ? weightx : NULL, &err,
168                                           distinct_callback, &ws[i]
169                                           );
170
171       for (; (c = casereader_read (rr)) != NULL; case_unref (c))
172         {
173           double sign = case_data (c, ws[i].sign)->f;
174           double rank = case_data_idx (c, weight ? 3 : 2)->f;
175           double w = 1.0;
176           if (weight)
177             w = case_data (c, weightx)->f;
178
179           if ( sign > 0 )
180             {
181               ws[i].positives.sum += rank * w;
182               ws[i].positives.n += w;
183             }
184           else if (sign < 0)
185             {
186               ws[i].negatives.sum += rank * w;
187               ws[i].negatives.n += w;
188             }
189           else
190             NOT_REACHED ();
191         }
192
193       casereader_destroy (rr);
194     }
195
196   casereader_destroy (input);
197
198   var_destroy (weightx);
199
200   show_ranks_box (ws, t2s, dict);
201   show_tests_box (ws, t2s, exact, timer);
202
203   for (i = 0 ; i < t2s->n_pairs; ++i )
204     {
205       var_destroy (ws[i].sign);
206       var_destroy (ws[i].absdiff);
207     }
208
209   free (ws);
210 }
211
212
213 \f
214
215 #include "gettext.h"
216 #define _(msgid) gettext (msgid)
217
218 static void
219 show_ranks_box (const struct wilcoxon_state *ws,
220                 const struct two_sample_test *t2s,
221                 const struct dictionary *dict)
222 {
223   size_t i;
224
225   const struct variable *wv = dict_get_weight (dict);
226   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : & F_8_0;
227
228   struct tab_table *table = tab_create (5, 1 + 4 * t2s->n_pairs);
229
230   tab_dim (table, tab_natural_dimensions, NULL, NULL);
231
232   tab_title (table, _("Ranks"));
233
234   tab_headers (table, 2, 0, 1, 0);
235
236   /* Vertical lines inside the box */
237   tab_box (table, 0, 0, -1, TAL_1,
238            1, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
239
240   /* Box around entire table */
241   tab_box (table, TAL_2, TAL_2, -1, -1,
242            0, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
243
244
245   tab_text (table,  2, 0,  TAB_CENTER, _("N"));
246   tab_text (table,  3, 0,  TAB_CENTER, _("Mean Rank"));
247   tab_text (table,  4, 0,  TAB_CENTER, _("Sum of Ranks"));
248
249
250   for (i = 0 ; i < t2s->n_pairs; ++i)
251     {
252       variable_pair *vp = &t2s->pairs[i];
253
254       struct string pair_name;
255       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
256       ds_put_cstr (&pair_name, " - ");
257       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
258
259       tab_text (table, 1, 1 + i * 4, TAB_LEFT, _("Negative Ranks"));
260       tab_text (table, 1, 2 + i * 4, TAB_LEFT, _("Positive Ranks"));
261       tab_text (table, 1, 3 + i * 4, TAB_LEFT, _("Ties"));
262       tab_text (table, 1, 4 + i * 4, TAB_LEFT, _("Total"));
263
264       tab_hline (table, TAL_1, 0, tab_nc (table) - 1, 1 + i * 4);
265
266
267       tab_text (table, 0, 1 + i * 4, TAB_LEFT, ds_cstr (&pair_name));
268       ds_destroy (&pair_name);
269
270
271       /* N */
272       tab_double (table, 2, 1 + i * 4, TAB_RIGHT, ws[i].negatives.n, wfmt);
273       tab_double (table, 2, 2 + i * 4, TAB_RIGHT, ws[i].positives.n, wfmt);
274       tab_double (table, 2, 3 + i * 4, TAB_RIGHT, ws[i].n_zeros, wfmt);
275
276       tab_double (table, 2, 4 + i * 4, TAB_RIGHT,
277                  ws[i].n_zeros + ws[i].positives.n + ws[i].negatives.n, wfmt);
278
279       /* Sums */
280       tab_double (table, 4, 1 + i * 4, TAB_RIGHT, ws[i].negatives.sum, NULL);
281       tab_double (table, 4, 2 + i * 4, TAB_RIGHT, ws[i].positives.sum, NULL);
282
283
284       /* Means */
285       tab_double (table, 3, 1 + i * 4, TAB_RIGHT,
286                  ws[i].negatives.sum / (double) ws[i].negatives.n, NULL);
287
288       tab_double (table, 3, 2 + i * 4, TAB_RIGHT,
289                  ws[i].positives.sum / (double) ws[i].positives.n, NULL);
290
291     }
292
293   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, 1);
294   tab_vline (table, TAL_2, 2, 0, tab_nr (table) - 1);
295
296
297   tab_submit (table);
298 }
299
300
301 static void
302 show_tests_box (const struct wilcoxon_state *ws,
303                 const struct two_sample_test *t2s,
304                 bool exact,
305                 double timer UNUSED
306                 )
307 {
308   size_t i;
309   struct tab_table *table = tab_create (1 + t2s->n_pairs, exact ? 5 : 3);
310
311   tab_dim (table, tab_natural_dimensions, NULL, NULL);
312
313   tab_title (table, _("Test Statistics"));
314
315   tab_headers (table, 1, 0, 1, 0);
316
317   /* Vertical lines inside the box */
318   tab_box (table, 0, 0, -1, TAL_1,
319            0, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
320
321   /* Box around entire 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_text (table,  0, 1,  TAB_LEFT, _("Z"));
327   tab_text (table,  0, 2,  TAB_LEFT, _("Asymp. Sig. (2-tailed)"));
328
329   if ( exact )
330     {
331       tab_text (table,  0, 3,  TAB_LEFT, _("Exact Sig. (2-tailed)"));
332       tab_text (table,  0, 4,  TAB_LEFT, _("Exact Sig. (1-tailed)"));
333
334 #if 0
335       tab_text (table,  0, 5,  TAB_LEFT, _("Point Probability"));
336 #endif
337     }
338
339   for (i = 0 ; i < t2s->n_pairs; ++i)
340     {
341       double z;
342       double n = ws[i].positives.n + ws[i].negatives.n;
343       variable_pair *vp = &t2s->pairs[i];
344
345       struct string pair_name;
346       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
347       ds_put_cstr (&pair_name, " - ");
348       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
349
350
351       tab_text (table, 1 + i, 0, TAB_CENTER, ds_cstr (&pair_name));
352       ds_destroy (&pair_name);
353
354       z = MIN (ws[i].positives.sum, ws[i].negatives.sum);
355       z -= n * (n + 1)/ 4.0;
356
357       z /= sqrt (n * (n + 1) * (2*n + 1)/24.0 - ws[i].tiebreaker / 48.0);
358
359       tab_double (table, 1 + i, 1, TAB_RIGHT, z, NULL);
360
361       tab_double (table, 1 + i, 2, TAB_RIGHT,
362                  2.0 * gsl_cdf_ugaussian_P (z),
363                  NULL);
364
365       if (exact)
366         {
367           double p = LevelOfSignificanceWXMPSR (ws[i].positives.sum, n);
368           if (p < 0)
369             {
370               msg (MW, ("Too many pairs to calculate exact significance."));
371             }
372           else
373             {
374               tab_double (table, 1 + i, 3, TAB_RIGHT, p, NULL);
375               tab_double (table, 1 + i, 4, TAB_RIGHT, p / 2.0, NULL);
376             }
377         }
378     }
379
380   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, 1);
381   tab_vline (table, TAL_2, 1, 0, tab_nr (table) - 1);
382
383
384   tab_submit (table);
385 }