pintos-os.org Git - pspp/blob - src/language/commands/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/commands/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_num (c, (*vp)[0]) - case_num (c, (*vp)[1]);
  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 caseproto *proto;
  91
  92   input =
  93     casereader_create_filter_weight (input, dict, &warn, NULL);
  94
  95   proto = caseproto_create ();
  96   proto = caseproto_add_width (proto, 0);
  97   proto = caseproto_add_width (proto, 0);
  98   if (weight != NULL)
  99     proto = caseproto_add_width (proto, 0);
 100
 101   for (i = 0 ; i < t2s->n_pairs; ++i)
 102     {
 103       struct casereader *r = casereader_clone (input);
 104       struct casewriter *writer;
 105       struct ccase *c;
 106       struct subcase ordering;
 107       variable_pair *vp = &t2s->pairs[i];
 108
 109       ws[i].dict = dict_create ("UTF-8");
 110       ws[i].sign = dict_create_var (ws[i].dict, "sign", 0);
 111       ws[i].absdiff = dict_create_var (ws[i].dict, "absdiff", 0);
 112       ws[i].weight = dict_create_var (ws[i].dict, "weight", 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_uninit (&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             *case_num_rw (output, ws[i].sign) = 1.0;
 129           else if (d < 0)
 130             *case_num_rw (output, ws[i].sign) = -1.0;
 131           else
 132             {
 133               double w = 1.0;
 134               if (weight)
 135                 w = case_num (c, weight);
 136
 137               /* Central point values should be dropped */
 138               ws[i].n_zeros += w;
 139               case_unref (output);
 140               continue;
 141             }
 142
 143           *case_num_rw (output, ws[i].absdiff) = fabs (d);
 144
 145           if (weight)
 146            *case_num_rw (output, ws[i].weight) = case_num (c, weight);
 147
 148           casewriter_write (writer, output);
 149         }
 150       casereader_destroy (r);
 151       ws[i].reader = casewriter_make_reader (writer);
 152     }
 153   caseproto_unref (proto);
 154
 155   for (i = 0 ; i < t2s->n_pairs; ++i)
 156     {
 157       struct casereader *rr ;
 158       struct ccase *c;
 159       enum rank_error err = 0;
 160
 161       rr = casereader_create_append_rank (ws[i].reader, ws[i].absdiff,
 162                                           weight ? ws[i].weight : NULL, &err,
 163                                           distinct_callback, &ws[i]
 164                                         );
 165
 166       for (; (c = casereader_read (rr)) != NULL; case_unref (c))
 167         {
 168           double sign = case_num (c, ws[i].sign);
 169           double rank = case_num_idx (c, weight ? 3 : 2);
 170           double w = weight ? case_num (c, ws[i].weight) : 1.0;
 171
 172           if (sign > 0)
 173             {
 174               ws[i].positives.sum += rank * w;
 175               ws[i].positives.n += w;
 176             }
 177           else if (sign < 0)
 178             {
 179               ws[i].negatives.sum += rank * w;
 180               ws[i].negatives.n += w;
 181             }
 182           else
 183             NOT_REACHED ();
 184         }
 185
 186       casereader_destroy (rr);
 187     }
 188
 189   casereader_destroy (input);
 190
 191   show_ranks_box (ws, t2s, dict);
 192   show_tests_box (ws, t2s, exact, timer);
 193
 194   for (i = 0 ; i < t2s->n_pairs; ++i)
 195     dict_unref (ws[i].dict);
 196
 197   free (ws);
 198 }
 199 \f
 200 static void
 201 put_row (struct pivot_table *table, int var_idx, int sign_idx,
 202          double n, double sum)
 203 {
 204   pivot_table_put3 (table, 0, sign_idx, var_idx, pivot_value_new_number (n));
 205   if (sum != SYSMIS)
 206     {
 207       pivot_table_put3 (table, 1, sign_idx, var_idx,
 208                         pivot_value_new_number (sum / n));
 209       pivot_table_put3 (table, 2, sign_idx, var_idx,
 210                         pivot_value_new_number (sum));
 211     }
 212 }
 213
 214 static int
 215 add_pair_leaf (struct pivot_dimension *dimension, variable_pair *pair)
 216 {
 217   char *label = xasprintf ("%s - %s", var_to_string ((*pair)[0]),
 218                            var_to_string ((*pair)[1]));
 219   return pivot_category_create_leaf (
 220     dimension->root,
 221     pivot_value_new_user_text_nocopy (label));
 222 }
 223
 224 static void
 225 show_ranks_box (const struct wilcoxon_state *ws,
 226                 const struct two_sample_test *t2s,
 227                 const struct dictionary *dict)
 228 {
 229   struct pivot_table *table = pivot_table_create (N_("Ranks"));
 230   pivot_table_set_weight_var (table, dict_get_weight (dict));
 231
 232   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 233                           N_("N"), PIVOT_RC_COUNT,
 234                           N_("Mean Rank"), PIVOT_RC_OTHER,
 235                           N_("Sum of Ranks"), PIVOT_RC_OTHER);
 236
 237   pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Sign"),
 238                           N_("Negative Ranks"), N_("Positive Ranks"),
 239                           N_("Ties"), N_("Total"));
 240
 241   struct pivot_dimension *pairs = pivot_dimension_create (
 242     table, PIVOT_AXIS_ROW, N_("Pairs"));
 243
 244   for (size_t i = 0 ; i < t2s->n_pairs; ++i)
 245     {
 246       variable_pair *vp = &t2s->pairs[i];
 247       int pair_idx = add_pair_leaf (pairs, vp);
 248
 249       const struct wilcoxon_state *w = &ws[i];
 250       put_row (table, pair_idx, 0, w->negatives.n, w->negatives.sum);
 251       put_row (table, pair_idx, 1, w->positives.n, w->positives.sum);
 252       put_row (table, pair_idx, 2, w->n_zeros, SYSMIS);
 253       put_row (table, pair_idx, 3,
 254                w->n_zeros + w->positives.n + w->negatives.n, SYSMIS);
 255     }
 256
 257   pivot_table_submit (table);
 258 }
 259
 260
 261 static void
 262 show_tests_box (const struct wilcoxon_state *ws,
 263                 const struct two_sample_test *t2s,
 264                 bool exact,
 265                 double timer UNUSED
 266                 )
 267 {
 268   struct pivot_table *table = pivot_table_create (N_("Test Statistics"));
 269
 270   struct pivot_dimension *statistics = pivot_dimension_create (
 271     table, PIVOT_AXIS_ROW, N_("Statistics"),
 272     N_("Z"), PIVOT_RC_OTHER,
 273     N_("Asymp. Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE);
 274   if (exact)
 275     pivot_category_create_leaves (
 276       statistics->root,
 277       N_("Exact Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE,
 278       N_("Exact Sig. (1-tailed)"), PIVOT_RC_SIGNIFICANCE);
 279
 280   struct pivot_dimension *pairs = pivot_dimension_create (
 281     table, PIVOT_AXIS_COLUMN, N_("Pairs"));
 282
 283   struct pivot_footnote *too_many_pairs = pivot_table_create_footnote (
 284     table, pivot_value_new_text (
 285       N_("Too many pairs to calculate exact significance")));
 286
 287   for (size_t i = 0 ; i < t2s->n_pairs; ++i)
 288     {
 289       variable_pair *vp = &t2s->pairs[i];
 290       int pair_idx = add_pair_leaf (pairs, vp);
 291
 292       double n = ws[i].positives.n + ws[i].negatives.n;
 293       double z = MIN (ws[i].positives.sum, ws[i].negatives.sum);
 294       z -= n * (n + 1)/ 4.0;
 295       z /= sqrt (n * (n + 1) * (2*n + 1)/24.0 - ws[i].tiebreaker / 48.0);
 296
 297       double entries[4];
 298       int n_entries = 0;
 299       entries[n_entries++] = z;
 300       entries[n_entries++] = 2.0 * gsl_cdf_ugaussian_P (z);
 301
 302       int footnote_idx = -1;
 303       if (exact)
 304         {
 305           double p = LevelOfSignificanceWXMPSR (ws[i].positives.sum, n);
 306           if (p < 0)
 307             {
 308               footnote_idx = n_entries;
 309               entries[n_entries++] = SYSMIS;
 310             }
 311           else
 312             {
 313               entries[n_entries++] = p;
 314               entries[n_entries++] = p / 2.0;
 315             }
 316         }
 317
 318       for (int j = 0; j < n_entries; j++)
 319         {
 320           struct pivot_value *value = pivot_value_new_number (entries[j]);
 321           if (j == footnote_idx)
 322             pivot_value_add_footnote (value, too_many_pairs);
 323           pivot_table_put2 (table, j, pair_idx, value);
 324         }
 325     }
 326
 327   pivot_table_submit (table);
 328 }