pintos-os.org Git - pspp/blob - src/language/stats/wilcoxon.c

   1 /* Pspp - a program for statistical analysis.
   2    Copyright (C) 2008, 2009, 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
  19 #include <config.h>
  20
  21 #include "language/stats/wilcoxon.h"
  22
  23 #include <gsl/gsl_cdf.h>
  24 #include <math.h>
  25
  26 #include "data/casereader.h"
  27 #include "data/casewriter.h"
  28 #include "data/dataset.h"
  29 #include "data/dictionary.h"
  30 #include "data/format.h"
  31 #include "data/subcase.h"
  32 #include "data/variable.h"
  33 #include "libpspp/assertion.h"
  34 #include "libpspp/message.h"
  35 #include "libpspp/misc.h"
  36 #include "math/sort.h"
  37 #include "math/wilcoxon-sig.h"
  38 #include "output/pivot-table.h"
  39
  40 #include "gl/minmax.h"
  41 #include "gl/xalloc.h"
  42
  43 #include "gettext.h"
  44 #define N_(msgid) msgid
  45 #define _(msgid) gettext (msgid)
  46
  47 static double
  48 append_difference (const struct ccase *c, casenumber n UNUSED, void *aux)
  49 {
  50   const variable_pair *vp = aux;
  51
  52   return case_data (c, (*vp)[0])->f - case_data (c, (*vp)[1])->f;
  53 }
  54
  55 static void show_ranks_box (const struct wilcoxon_state *,
  56                             const struct two_sample_test *,
  57                             const struct dictionary *);
  58
  59 static void show_tests_box (const struct wilcoxon_state *,
  60                             const struct two_sample_test *,
  61                             bool exact, double timer);
  62
  63
  64
  65 static void
  66 distinct_callback (double v UNUSED, casenumber n, double w UNUSED, void *aux)
  67 {
  68   struct wilcoxon_state *ws = aux;
  69
  70   ws->tiebreaker += pow3 (n) - n;
  71 }
  72
  73 #define WEIGHT_IDX 2
  74
  75 void
  76 wilcoxon_execute (const struct dataset *ds,
  77                   struct casereader *input,
  78                   enum mv_class exclude,
  79                   const struct npar_test *test,
  80                   bool exact,
  81                   double timer)
  82 {
  83   int i;
  84   bool warn = true;
  85   const struct dictionary *dict = dataset_dict (ds);
  86   const struct two_sample_test *t2s = UP_CAST (test, const struct two_sample_test, parent);
  87
  88   struct wilcoxon_state *ws = XCALLOC (t2s->n_pairs,  struct wilcoxon_state);
  89   const struct variable *weight = dict_get_weight (dict);
  90   struct variable *weightx = dict_create_internal_var (WEIGHT_IDX, 0);
  91   struct caseproto *proto;
  92
  93   input =
  94     casereader_create_filter_weight (input, dict, &warn, NULL);
  95
  96   proto = caseproto_create ();
  97   proto = caseproto_add_width (proto, 0);
  98   proto = caseproto_add_width (proto, 0);
  99   if (weight != NULL)
 100     proto = caseproto_add_width (proto, 0);
 101
 102   for (i = 0 ; i < t2s->n_pairs; ++i)
 103     {
 104       struct casereader *r = casereader_clone (input);
 105       struct casewriter *writer;
 106       struct ccase *c;
 107       struct subcase ordering;
 108       variable_pair *vp = &t2s->pairs[i];
 109
 110       ws[i].sign = dict_create_internal_var (0, 0);
 111       ws[i].absdiff = dict_create_internal_var (1, 0);
 112
 113       r = casereader_create_filter_missing (r, *vp, 2,
 114                                             exclude,
 115                                             NULL, NULL);
 116
 117       subcase_init_var (&ordering, ws[i].absdiff, SC_ASCEND);
 118       writer = sort_create_writer (&ordering, proto);
 119       subcase_destroy (&ordering);
 120
 121       for (; (c = casereader_read (r)) != NULL; case_unref (c))
 122         {
 123           struct ccase *output = case_create (proto);
 124           double d = append_difference (c, 0, vp);
 125
 126           if (d > 0)
 127             {
 128               case_data_rw (output, ws[i].sign)->f = 1.0;
 129
 130             }
 131           else if (d < 0)
 132             {
 133               case_data_rw (output, ws[i].sign)->f = -1.0;
 134             }
 135           else
 136             {
 137               double w = 1.0;
 138               if (weight)
 139                 w = case_data (c, weight)->f;
 140
 141               /* Central point values should be dropped */
 142               ws[i].n_zeros += w;
 143               case_unref (output);
 144               continue;
 145             }
 146
 147           case_data_rw (output, ws[i].absdiff)->f = fabs (d);
 148
 149           if (weight)
 150            case_data_rw (output, weightx)->f = case_data (c, weight)->f;
 151
 152           casewriter_write (writer, output);
 153         }
 154       casereader_destroy (r);
 155       ws[i].reader = casewriter_make_reader (writer);
 156     }
 157   caseproto_unref (proto);
 158
 159   for (i = 0 ; i < t2s->n_pairs; ++i)
 160     {
 161       struct casereader *rr ;
 162       struct ccase *c;
 163       enum rank_error err = 0;
 164
 165       rr = casereader_create_append_rank (ws[i].reader, ws[i].absdiff,
 166                                           weight ? weightx : NULL, &err,
 167                                           distinct_callback, &ws[i]
 168                                         );
 169
 170       for (; (c = casereader_read (rr)) != NULL; case_unref (c))
 171         {
 172           double sign = case_data (c, ws[i].sign)->f;
 173           double rank = case_data_idx (c, weight ? 3 : 2)->f;
 174           double w = 1.0;
 175           if (weight)
 176             w = case_data (c, weightx)->f;
 177
 178           if (sign > 0)
 179             {
 180               ws[i].positives.sum += rank * w;
 181               ws[i].positives.n += w;
 182             }
 183           else if (sign < 0)
 184             {
 185               ws[i].negatives.sum += rank * w;
 186               ws[i].negatives.n += w;
 187             }
 188           else
 189             NOT_REACHED ();
 190         }
 191
 192       casereader_destroy (rr);
 193     }
 194
 195   casereader_destroy (input);
 196
 197   dict_destroy_internal_var (weightx);
 198
 199   show_ranks_box (ws, t2s, dict);
 200   show_tests_box (ws, t2s, exact, timer);
 201
 202   for (i = 0 ; i < t2s->n_pairs; ++i)
 203     {
 204       dict_destroy_internal_var (ws[i].sign);
 205       dict_destroy_internal_var (ws[i].absdiff);
 206     }
 207
 208   free (ws);
 209 }
 210 \f
 211 static void
 212 put_row (struct pivot_table *table, int var_idx, int sign_idx,
 213          double n, double sum)
 214 {
 215   pivot_table_put3 (table, 0, sign_idx, var_idx, pivot_value_new_number (n));
 216   if (sum != SYSMIS)
 217     {
 218       pivot_table_put3 (table, 1, sign_idx, var_idx,
 219                         pivot_value_new_number (sum / n));
 220       pivot_table_put3 (table, 2, sign_idx, var_idx,
 221                         pivot_value_new_number (sum));
 222     }
 223 }
 224
 225 static int
 226 add_pair_leaf (struct pivot_dimension *dimension, variable_pair *pair)
 227 {
 228   char *label = xasprintf ("%s - %s", var_to_string ((*pair)[0]),
 229                            var_to_string ((*pair)[1]));
 230   return pivot_category_create_leaf (
 231     dimension->root,
 232     pivot_value_new_user_text_nocopy (label));
 233 }
 234
 235 static void
 236 show_ranks_box (const struct wilcoxon_state *ws,
 237                 const struct two_sample_test *t2s,
 238                 const struct dictionary *dict)
 239 {
 240   struct pivot_table *table = pivot_table_create (N_("Ranks"));
 241   pivot_table_set_weight_var (table, dict_get_weight (dict));
 242
 243   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 244                           N_("N"), PIVOT_RC_COUNT,
 245                           N_("Mean Rank"), PIVOT_RC_OTHER,
 246                           N_("Sum of Ranks"), PIVOT_RC_OTHER);
 247
 248   pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Sign"),
 249                           N_("Negative Ranks"), N_("Positive Ranks"),
 250                           N_("Ties"), N_("Total"));
 251
 252   struct pivot_dimension *pairs = pivot_dimension_create (
 253     table, PIVOT_AXIS_ROW, N_("Pairs"));
 254
 255   for (size_t i = 0 ; i < t2s->n_pairs; ++i)
 256     {
 257       variable_pair *vp = &t2s->pairs[i];
 258       int pair_idx = add_pair_leaf (pairs, vp);
 259
 260       const struct wilcoxon_state *w = &ws[i];
 261       put_row (table, pair_idx, 0, w->negatives.n, w->negatives.sum);
 262       put_row (table, pair_idx, 1, w->positives.n, w->positives.sum);
 263       put_row (table, pair_idx, 2, w->n_zeros, SYSMIS);
 264       put_row (table, pair_idx, 3,
 265                w->n_zeros + w->positives.n + w->negatives.n, SYSMIS);
 266     }
 267
 268   pivot_table_submit (table);
 269 }
 270
 271
 272 static void
 273 show_tests_box (const struct wilcoxon_state *ws,
 274                 const struct two_sample_test *t2s,
 275                 bool exact,
 276                 double timer UNUSED
 277                 )
 278 {
 279   struct pivot_table *table = pivot_table_create (N_("Test Statistics"));
 280
 281   struct pivot_dimension *statistics = pivot_dimension_create (
 282     table, PIVOT_AXIS_ROW, N_("Statistics"),
 283     N_("Z"), PIVOT_RC_OTHER,
 284     N_("Asymp. Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE);
 285   if (exact)
 286     pivot_category_create_leaves (
 287       statistics->root,
 288       N_("Exact Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE,
 289       N_("Exact Sig. (1-tailed)"), PIVOT_RC_SIGNIFICANCE);
 290
 291   struct pivot_dimension *pairs = pivot_dimension_create (
 292     table, PIVOT_AXIS_COLUMN, N_("Pairs"));
 293
 294   struct pivot_footnote *too_many_pairs = pivot_table_create_footnote (
 295     table, pivot_value_new_text (
 296       N_("Too many pairs to calculate exact significance")));
 297
 298   for (size_t i = 0 ; i < t2s->n_pairs; ++i)
 299     {
 300       variable_pair *vp = &t2s->pairs[i];
 301       int pair_idx = add_pair_leaf (pairs, vp);
 302
 303       double n = ws[i].positives.n + ws[i].negatives.n;
 304       double z = MIN (ws[i].positives.sum, ws[i].negatives.sum);
 305       z -= n * (n + 1)/ 4.0;
 306       z /= sqrt (n * (n + 1) * (2*n + 1)/24.0 - ws[i].tiebreaker / 48.0);
 307
 308       double entries[4];
 309       int n_entries = 0;
 310       entries[n_entries++] = z;
 311       entries[n_entries++] = 2.0 * gsl_cdf_ugaussian_P (z);
 312
 313       int footnote_idx = -1;
 314       if (exact)
 315         {
 316           double p = LevelOfSignificanceWXMPSR (ws[i].positives.sum, n);
 317           if (p < 0)
 318             {
 319               footnote_idx = n_entries;
 320               entries[n_entries++] = SYSMIS;
 321             }
 322           else
 323             {
 324               entries[n_entries++] = p;
 325               entries[n_entries++] = p / 2.0;
 326             }
 327         }
 328
 329       for (int j = 0; j < n_entries; j++)
 330         {
 331           struct pivot_value *value = pivot_value_new_number (entries[j]);
 332           if (j == footnote_idx)
 333             pivot_value_add_footnote (value, too_many_pairs);
 334           pivot_table_put2 (table, j, pair_idx, value);
 335         }
 336     }
 337
 338   pivot_table_submit (table);
 339 }