src/language/stats/friedman.c

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