pintos-os.org Git - pspp/blob - 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/pivot-table.h"
  33
  34 #include "gettext.h"
  35 #define N_(msgid) msgid
  36 #define _(msgid) gettext (msgid)
  37
  38
  39 struct friedman
  40 {
  41   double *rank_sum;
  42   double cc;
  43   double chi_sq;
  44   double w;
  45   const struct dictionary *dict;
  46 };
  47
  48 static void show_ranks_box (const struct one_sample_test *ost,
  49                             const struct friedman *fr);
  50
  51 static void show_sig_box (const struct one_sample_test *ost,
  52                           const struct friedman *fr);
  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 *ost = UP_CAST (test, struct one_sample_test, parent);
 100   struct friedman_test *ft = UP_CAST (ost, struct friedman_test, parent);
 101   bool warn = true;
 102
 103   double sigma_t = 0.0;
 104   struct datum *row = xcalloc (ost->n_vars, sizeof *row);
 105   double rsq;
 106   struct friedman fr;
 107   fr.rank_sum = xcalloc (ost->n_vars, sizeof *fr.rank_sum);
 108   fr.cc = 0.0;
 109   fr.dict = dict;
 110   for (v = 0; v < ost->n_vars; ++v)
 111     {
 112       row[v].posn = v;
 113       fr.rank_sum[v] = 0.0;
 114     }
 115
 116   input = casereader_create_filter_weight (input, dict, &warn, NULL);
 117   input = casereader_create_filter_missing (input,
 118                                             ost->vars, ost->n_vars,
 119                                             exclude, 0, 0);
 120
 121   for (; (c = casereader_read (input)); case_unref (c))
 122     {
 123       double prev_x = SYSMIS;
 124       int run_length = 0;
 125
 126       const double w = weight ? case_data (c, weight)->f: 1.0;
 127
 128       fr.cc += w;
 129
 130       for (v = 0; v < ost->n_vars; ++v)
 131         {
 132           const struct variable *var = ost->vars[v];
 133           const union value *val = case_data (c, var);
 134           row[v].x = val->f;
 135         }
 136
 137       qsort (row, ost->n_vars, sizeof *row, cmp_x);
 138       for (v = 0; v < ost->n_vars; ++v)
 139         {
 140           double x = row[v].x;
 141           /* Replace value by the Rank */
 142           if ( prev_x == x)
 143             {
 144               /* Deal with ties */
 145               int i;
 146               run_length++;
 147               for (i = v - run_length; i < v; ++i)
 148                 {
 149                   row[i].x *= run_length ;
 150                   row[i].x += v + 1;
 151                   row[i].x /= run_length + 1;
 152                 }
 153               row[v].x = row[v-1].x;
 154             }
 155           else
 156             {
 157               row[v].x = v + 1;
 158               if ( run_length > 0)
 159                 {
 160                   double t = run_length + 1;
 161                   sigma_t += w * (pow3 (t) - t);
 162                 }
 163               run_length = 0;
 164             }
 165           prev_x = x;
 166         }
 167       if ( run_length > 0)
 168         {
 169           double t = run_length + 1;
 170           sigma_t += w * (pow3 (t) - t );
 171         }
 172
 173       qsort (row, ost->n_vars, sizeof *row, cmp_posn);
 174
 175       for (v = 0; v < ost->n_vars; ++v)
 176         fr.rank_sum[v] += row[v].x * w;
 177     }
 178   casereader_destroy (input);
 179   free (row);
 180
 181
 182   for (v = 0; v < ost->n_vars; ++v)
 183     {
 184       numerator += pow2 (fr.rank_sum[v]);
 185     }
 186
 187   rsq = numerator;
 188
 189   numerator *= 12.0 / (fr.cc * ost->n_vars * ( ost->n_vars + 1));
 190   numerator -= 3 * fr.cc * ( ost->n_vars + 1);
 191
 192   denominator = 1 - sigma_t / ( fr.cc * ost->n_vars * ( pow2 (ost->n_vars) - 1));
 193
 194   fr.chi_sq = numerator / denominator;
 195
 196   if ( ft->kendalls_w)
 197     {
 198       fr.w = 12 * rsq ;
 199       fr.w -= 3 * pow2 (fr.cc) *
 200         ost->n_vars * pow2 (ost->n_vars + 1);
 201
 202       fr.w /= pow2 (fr.cc) * (pow3 (ost->n_vars) - ost->n_vars)
 203         - fr.cc * sigma_t;
 204     }
 205   else
 206     fr.w = SYSMIS;
 207
 208   show_ranks_box (ost, &fr);
 209   show_sig_box (ost, &fr);
 210
 211   free (fr.rank_sum);
 212 }
 213
 214 \f
 215
 216
 217 static void
 218 show_ranks_box (const struct one_sample_test *ost, const struct friedman *fr)
 219 {
 220   struct pivot_table *table = pivot_table_create (N_("Ranks"));
 221
 222   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Mean Rank"),
 223                           N_("Mean Rank"), PIVOT_RC_OTHER);
 224
 225   struct pivot_dimension *variables = pivot_dimension_create (
 226     table, PIVOT_AXIS_ROW, N_("Variable"));
 227
 228   for (size_t i = 0 ; i < ost->n_vars ; ++i)
 229     {
 230       int row = pivot_category_create_leaf (
 231         variables->root, pivot_value_new_variable (ost->vars[i]));
 232
 233       pivot_table_put2 (table, 0, row,
 234                         pivot_value_new_number (fr->rank_sum[i] / fr->cc));
 235     }
 236
 237   pivot_table_submit (table);
 238 }
 239
 240
 241 static void
 242 show_sig_box (const struct one_sample_test *ost, const struct friedman *fr)
 243 {
 244   const struct friedman_test *ft = UP_CAST (ost, const struct friedman_test, parent);
 245
 246   struct pivot_table *table = pivot_table_create (N_("Test Statistics"));
 247   pivot_table_set_weight_var (table, dict_get_weight (fr->dict));
 248
 249   struct pivot_dimension *statistics = pivot_dimension_create (
 250     table, PIVOT_AXIS_ROW, N_("Statistics"),
 251     N_("N"), PIVOT_RC_COUNT);
 252   if (ft->kendalls_w)
 253     pivot_category_create_leaves (statistics->root, N_("Kendall's W"),
 254                                   PIVOT_RC_OTHER);
 255   pivot_category_create_leaves (statistics->root,
 256                                 N_("Chi-Square"), PIVOT_RC_OTHER,
 257                                 N_("df"), PIVOT_RC_INTEGER,
 258                                 N_("Asymp. Sig."), PIVOT_RC_SIGNIFICANCE);
 259
 260   double entries[5];
 261   int n = 0;
 262
 263   entries[n++] = fr->cc;
 264   if (ft->kendalls_w)
 265     entries[n++] = fr->w;
 266   entries[n++] = fr->chi_sq;
 267   entries[n++] = ost->n_vars - 1;
 268   entries[n++] = gsl_cdf_chisq_Q (fr->chi_sq, ost->n_vars - 1);
 269   assert (n <= sizeof entries / sizeof *entries);
 270
 271   for (size_t i = 0; i < n; i++)
 272     pivot_table_put1 (table, i, pivot_value_new_number (entries[i]));
 273
 274   pivot_table_submit (table);
 275 }