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