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/case-ordering.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 case_ordering *ordering = case_ordering_create ();
  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       case_ordering_add_var (ordering, ws[i].absdiff, SRT_ASCEND);
 101
 102
 103       r = casereader_create_filter_missing (r, *vp, 2,
 104                                             exclude,
 105                                             NULL, NULL);
 106
 107       writer = sort_create_writer (ordering, reader_width);
 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_restorer = 0;
 406   timeout_action.sa_handler = give_up_callback;
 407
 408   if ( 0 == sigsetjmp (env, 1))
 409     {
 410       sigaction (SIGALRM, &timeout_action, &old_action);
 411       alarm (timer * 60.0);
 412
 413       p = LevelOfSignificanceWXMPSR (w, n);
 414     }
 415
 416   sigaction (SIGALRM, &old_action, NULL);
 417
 418   return p;
 419 }