src/language/stats/friedman.c

   1 /* PSPP - a program for statistical analysis. -*-c-*-
   2    Copyright (C) 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 #include <config.h>
  19
  20 #include "language/stats/friedman.h"
  21
  22 #include <gsl/gsl_cdf.h>
  23 #include <math.h>
  24
  25 #include "data/casereader.h"
  26 #include "data/dataset.h"
  27 #include "data/dictionary.h"
  28 #include "data/format.h"
  29 #include "data/variable.h"
  30 #include "libpspp/message.h"
  31 #include "libpspp/misc.h"
  32 #include "output/tab.h"
  33
  34 #include "gettext.h"
  35 #define _(msgid) gettext (msgid)
  36
  37
  38 struct friedman
  39 {
  40   double *rank_sum;
  41   double cc;
  42   double chi_sq;
  43   double w;
  44   const struct dictionary *dict;
  45 };
  46
  47 static void show_ranks_box (const struct one_sample_test *ost,
  48                             const struct friedman *fr);
  49
  50 static void show_sig_box (const struct one_sample_test *ost,
  51                           const struct friedman *fr);
  52
  53 struct datum
  54 {
  55   long posn;
  56   double x;
  57 };
  58
  59 static int
  60 cmp_x (const void *a_, const void *b_)
  61 {
  62   const struct datum *a = a_;
  63   const struct datum *b = b_;
  64
  65   if (a->x < b->x)
  66     return -1;
  67
  68   return (a->x > b->x);
  69 }
  70
  71 static int
  72 cmp_posn (const void *a_, const void *b_)
  73 {
  74   const struct datum *a = a_;
  75   const struct datum *b = b_;
  76
  77   if (a->posn < b->posn)
  78     return -1;
  79
  80   return (a->posn > b->posn);
  81 }
  82
  83 void
  84 friedman_execute (const struct dataset *ds,
  85                   struct casereader *input,
  86                   enum mv_class exclude,
  87                   const struct npar_test *test,
  88                   bool exact UNUSED,
  89                   double timer UNUSED)
  90 {
  91   double numerator = 0.0;
  92   double denominator = 0.0;
  93   int v;
  94   struct ccase *c;
  95   const struct dictionary *dict = dataset_dict (ds);
  96   const struct variable *weight = dict_get_weight (dict);
  97
  98   struct one_sample_test *ost = UP_CAST (test, struct one_sample_test, parent);
  99   struct friedman_test *ft = UP_CAST (ost, struct friedman_test, parent);
 100   bool warn = true;
 101
 102   double sigma_t = 0.0;
 103   struct datum *row = xcalloc (ost->n_vars, sizeof *row);
 104   double rsq;
 105   struct friedman fr;
 106   fr.rank_sum = xcalloc (ost->n_vars, sizeof *fr.rank_sum);
 107   fr.cc = 0.0;
 108   fr.dict = dict;
 109   for (v = 0; v < ost->n_vars; ++v)
 110     {
 111       row[v].posn = v;
 112       fr.rank_sum[v] = 0.0;
 113     }
 114
 115   input = casereader_create_filter_weight (input, dict, &warn, NULL);
 116   input = casereader_create_filter_missing (input,
 117                                             ost->vars, ost->n_vars,
 118                                             exclude, 0, 0);
 119
 120   for (; (c = casereader_read (input)); case_unref (c))
 121     {
 122       double prev_x = SYSMIS;
 123       int run_length = 0;
 124
 125       const double w = weight ? case_data (c, weight)->f: 1.0;
 126
 127       fr.cc += w;
 128
 129       for (v = 0; v < ost->n_vars; ++v)
 130         {
 131           const struct variable *var = ost->vars[v];
 132           const union value *val = case_data (c, var);
 133           row[v].x = val->f;
 134         }
 135
 136       qsort (row, ost->n_vars, sizeof *row, cmp_x);
 137       for (v = 0; v < ost->n_vars; ++v)
 138         {
 139           double x = row[v].x;
 140           /* Replace value by the Rank */
 141           if ( prev_x == x)
 142             {
 143               /* Deal with ties */
 144               int i;
 145               run_length++;
 146               for (i = v - run_length; i < v; ++i)
 147                 {
 148                   row[i].x *= run_length ;
 149                   row[i].x += v + 1;
 150                   row[i].x /= run_length + 1;
 151                 }
 152               row[v].x = row[v-1].x;
 153             }
 154           else
 155             {
 156               row[v].x = v + 1;
 157               if ( run_length > 0)
 158                 {
 159                   double t = run_length + 1;
 160                   sigma_t += w * (pow3 (t) - t);
 161                 }
 162               run_length = 0;
 163             }
 164           prev_x = x;
 165         }
 166       if ( run_length > 0)
 167         {
 168           double t = run_length + 1;
 169           sigma_t += w * (pow3 (t) - t );
 170         }
 171
 172       qsort (row, ost->n_vars, sizeof *row, cmp_posn);
 173
 174       for (v = 0; v < ost->n_vars; ++v)
 175         fr.rank_sum[v] += row[v].x * w;
 176     }
 177   casereader_destroy (input);
 178   free (row);
 179
 180
 181   for (v = 0; v < ost->n_vars; ++v)
 182     {
 183       numerator += pow2 (fr.rank_sum[v]);
 184     }
 185
 186   rsq = numerator;
 187
 188   numerator *= 12.0 / (fr.cc * ost->n_vars * ( ost->n_vars + 1));
 189   numerator -= 3 * fr.cc * ( ost->n_vars + 1);
 190
 191   denominator = 1 - sigma_t / ( fr.cc * ost->n_vars * ( pow2 (ost->n_vars) - 1));
 192
 193   fr.chi_sq = numerator / denominator;
 194
 195   if ( ft->kendalls_w)
 196     {
 197       fr.w = 12 * rsq ;
 198       fr.w -= 3 * pow2 (fr.cc) *
 199         ost->n_vars * pow2 (ost->n_vars + 1);
 200
 201       fr.w /= pow2 (fr.cc) * (pow3 (ost->n_vars) - ost->n_vars)
 202         - fr.cc * sigma_t;
 203     }
 204   else
 205     fr.w = SYSMIS;
 206
 207   show_ranks_box (ost, &fr);
 208   show_sig_box (ost, &fr);
 209
 210   free (fr.rank_sum);
 211 }
 212
 213 \f
 214
 215
 216 static void
 217 show_ranks_box (const struct one_sample_test *ost, const struct friedman *fr)
 218 {
 219   int i;
 220   const int row_headers = 1;
 221   const int column_headers = 1;
 222   struct tab_table *table =
 223     tab_create (row_headers + 1, column_headers + ost->n_vars);
 224
 225   tab_headers (table, row_headers, 0, column_headers, 0);
 226
 227   tab_title (table, _("Ranks"));
 228
 229   /* Vertical lines inside the box */
 230   tab_box (table, 1, 0, -1, TAL_1,
 231            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
 232
 233   /* Box around the table */
 234   tab_box (table, TAL_2, TAL_2, -1, -1,
 235            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 236
 237
 238   tab_text (table, 1, 0, 0, _("Mean Rank"));
 239
 240   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, column_headers);
 241   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 242
 243   for (i = 0 ; i < ost->n_vars ; ++i)
 244     {
 245       tab_text (table, 0, row_headers + i,
 246                 TAB_LEFT, var_to_string (ost->vars[i]));
 247
 248       tab_double (table, 1, row_headers + i,
 249                   0, fr->rank_sum[i] / fr->cc, 0);
 250     }
 251
 252   tab_submit (table);
 253 }
 254
 255
 256 static void
 257 show_sig_box (const struct one_sample_test *ost, const struct friedman *fr)
 258 {
 259   const struct friedman_test *ft = UP_CAST (ost, const struct friedman_test, parent);
 260
 261   int row = 0;
 262   const struct variable *weight = dict_get_weight (fr->dict);
 263   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 264
 265   const int row_headers = 1;
 266   const int column_headers = 0;
 267   struct tab_table *table =
 268     tab_create (row_headers + 1, column_headers + (ft->kendalls_w ? 5 : 4));
 269
 270   tab_headers (table, row_headers, 0, column_headers, 0);
 271
 272   tab_title (table, _("Test Statistics"));
 273
 274   tab_text (table,  0, column_headers + row++,
 275             TAT_TITLE | TAB_LEFT , _("N"));
 276
 277   if ( ft->kendalls_w)
 278     tab_text (table,  0, column_headers + row++,
 279               TAT_TITLE | TAB_LEFT , _("Kendall's W"));
 280
 281   tab_text (table,  0, column_headers + row++,
 282             TAT_TITLE | TAB_LEFT , _("Chi-Square"));
 283
 284   tab_text (table,  0, column_headers + row++,
 285             TAT_TITLE | TAB_LEFT, _("df"));
 286
 287   tab_text (table,  0, column_headers + row++,
 288             TAT_TITLE | TAB_LEFT, _("Asymp. Sig."));
 289
 290   /* Box around the table */
 291   tab_box (table, TAL_2, TAL_2, -1, -1,
 292            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 293
 294
 295   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 296   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 297
 298   row = 0;
 299   tab_double (table, 1, column_headers + row++,
 300               0, fr->cc, wfmt);
 301
 302   if (ft->kendalls_w)
 303     tab_double (table, 1, column_headers + row++,
 304                 0, fr->w, 0);
 305
 306   tab_double (table, 1, column_headers + row++,
 307               0, fr->chi_sq, 0);
 308
 309   tab_double (table, 1, column_headers + row++,
 310               0, ost->n_vars - 1, &F_8_0);
 311
 312   tab_double (table, 1, column_headers + row++,
 313               0, gsl_cdf_chisq_Q (fr->chi_sq, ost->n_vars - 1),
 314               0);
 315
 316   tab_submit (table);
 317 }