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

   1 /* Pspp - a program for statistical analysis.
   2    Copyright (C) 2010, 2011, 2022 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 "kruskal-wallis.h"
  21
  22 #include <gsl/gsl_cdf.h>
  23 #include <math.h>
  24
  25 #include "data/casereader.h"
  26 #include "data/casewriter.h"
  27 #include "data/dataset.h"
  28 #include "data/dictionary.h"
  29 #include "data/format.h"
  30 #include "data/subcase.h"
  31 #include "data/variable.h"
  32 #include "libpspp/assertion.h"
  33 #include "libpspp/hmap.h"
  34 #include "libpspp/bt.h"
  35 #include "libpspp/message.h"
  36 #include "libpspp/misc.h"
  37 #include "math/sort.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 /* Returns true iff the independent variable lies between nst->val1 and  nst->val2 */
  48 static bool
  49 include_func (const struct ccase *c, void *aux)
  50 {
  51   const struct n_sample_test *nst = aux;
  52
  53   const union value *smaller = 0;
  54   const union value *larger = 0;
  55   int x = value_compare_3way (&nst->val1, &nst->val2, var_get_width (nst->indep_var));
  56    if (x < 0)
  57     {
  58       smaller = &nst->val1;
  59       larger = &nst->val2;
  60     }
  61   else
  62     {
  63       smaller = &nst->val2;
  64       larger = &nst->val1;
  65     }
  66
  67   if (0 < value_compare_3way (smaller, case_data (c, nst->indep_var),
  68                               var_get_width (nst->indep_var)))
  69     return false;
  70
  71   if (0 > value_compare_3way (larger, case_data (c, nst->indep_var),
  72                               var_get_width (nst->indep_var)))
  73     return false;
  74
  75   return true;
  76 }
  77
  78
  79 struct rank_entry
  80 {
  81   struct hmap_node node;
  82   struct bt_node btn;
  83   union value group;
  84
  85   double sum_of_ranks;
  86   double n;
  87 };
  88
  89
  90 static int
  91 compare_rank_entries_3way (const struct bt_node *a,
  92                            const struct bt_node *b,
  93                            const void *aux)
  94 {
  95   const struct variable *var = aux;
  96   const struct rank_entry *rea = BT_DATA (a, struct rank_entry, btn);
  97   const struct rank_entry *reb = BT_DATA (b, struct rank_entry, btn);
  98
  99   return value_compare_3way (&rea->group, &reb->group, var_get_width (var));
 100 }
 101
 102
 103 /* Return the entry with the key GROUP or null if there is no such entry */
 104 static struct rank_entry *
 105 find_rank_entry (const struct hmap *map, const union value *group, size_t width)
 106 {
 107   struct rank_entry *re = NULL;
 108   size_t hash  = value_hash (group, width, 0);
 109
 110   HMAP_FOR_EACH_WITH_HASH (re, struct rank_entry, node, hash, map)
 111     {
 112       if (0 == value_compare_3way (group, &re->group, width))
 113         return re;
 114     }
 115
 116   return re;
 117 }
 118
 119 /* Calculates the adjustment necessary for tie compensation */
 120 static void
 121 distinct_callback (double v UNUSED, casenumber t, double w UNUSED, void *aux)
 122 {
 123   double *tiebreaker = aux;
 124
 125   *tiebreaker += pow3 (t) - t;
 126 }
 127
 128
 129 struct kw
 130 {
 131   struct hmap map;
 132   double h;
 133 };
 134
 135 static void show_ranks_box (const struct n_sample_test *, const struct kw *);
 136 static void show_sig_box (const struct n_sample_test *, const struct kw *);
 137
 138 void
 139 kruskal_wallis_execute (const struct dataset *ds,
 140                         struct casereader *input,
 141                         enum mv_class exclude,
 142                         const struct npar_test *test,
 143                         bool exact UNUSED,
 144                         double timer UNUSED)
 145 {
 146   int i;
 147   struct ccase *c;
 148   bool warn = true;
 149   const struct dictionary *dict = dataset_dict (ds);
 150   const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test, parent);
 151   const struct caseproto *proto ;
 152   size_t rank_idx ;
 153
 154   int total_n_groups = 0.0;
 155
 156   struct kw *kw = XCALLOC (nst->n_vars,  struct kw);
 157
 158   /* If the independent variable is missing, then we ignore the case */
 159   input = casereader_create_filter_missing (input,
 160                                             &nst->indep_var, 1,
 161                                             exclude,
 162                                             NULL, NULL);
 163
 164   input = casereader_create_filter_weight (input, dict, &warn, NULL);
 165
 166   /* Remove all those cases which are outside the range (val1, val2) */
 167   input = casereader_create_filter_func (input, include_func, NULL,
 168         CONST_CAST (struct n_sample_test *, nst), NULL);
 169
 170   proto = casereader_get_proto (input);
 171   rank_idx = caseproto_get_n_widths (proto);
 172
 173   /* Rank cases by the v value */
 174   for (i = 0; i < nst->n_vars; ++i)
 175     {
 176       double tiebreaker = 0.0;
 177       bool warn = true;
 178       enum rank_error rerr = 0;
 179       struct casereader *rr;
 180       struct casereader *r = casereader_clone (input);
 181
 182       r = sort_execute_1var (r, nst->vars[i]);
 183
 184       /* Ignore missings in the test variable */
 185       r = casereader_create_filter_missing (r, &nst->vars[i], 1,
 186                                             exclude,
 187                                             NULL, NULL);
 188
 189       rr = casereader_create_append_rank (r,
 190                                           nst->vars[i],
 191                                           dict_get_weight (dict),
 192                                           &rerr,
 193                                           distinct_callback, &tiebreaker);
 194
 195       hmap_init (&kw[i].map);
 196       for (; (c = casereader_read (rr)); case_unref (c))
 197         {
 198           const union value *group = case_data (c, nst->indep_var);
 199           const size_t group_var_width = var_get_width (nst->indep_var);
 200           struct rank_entry *rank = find_rank_entry (&kw[i].map, group, group_var_width);
 201
 202           if (NULL == rank)
 203             {
 204               rank = xzalloc (sizeof *rank);
 205               value_clone (&rank->group, group, group_var_width);
 206
 207               hmap_insert (&kw[i].map, &rank->node,
 208                            value_hash (&rank->group, group_var_width, 0));
 209             }
 210
 211           rank->sum_of_ranks += case_num_idx (c, rank_idx);
 212           rank->n += dict_get_case_weight (dict, c, &warn);
 213
 214           /* If this assertion fires, then either the data wasn't sorted or some other
 215              problem occurred */
 216           assert (rerr == 0);
 217         }
 218
 219       casereader_destroy (rr);
 220
 221       /* Calculate the value of h */
 222       {
 223         struct rank_entry *mre;
 224         double n = 0.0;
 225
 226         HMAP_FOR_EACH (mre, struct rank_entry, node, &kw[i].map)
 227           {
 228             kw[i].h += pow2 (mre->sum_of_ranks) / mre->n;
 229             n += mre->n;
 230
 231             total_n_groups ++;
 232           }
 233         kw[i].h *= 12 / (n * (n + 1));
 234         kw[i].h -= 3 * (n + 1) ;
 235
 236         kw[i].h /= 1 - tiebreaker/ (pow3 (n) - n);
 237       }
 238     }
 239
 240   casereader_destroy (input);
 241
 242   show_ranks_box (nst, kw);
 243   show_sig_box (nst, kw);
 244
 245   /* Cleanup allocated memory */
 246   for (i = 0 ; i < nst->n_vars; ++i)
 247     {
 248       struct rank_entry *mre, *next;
 249       HMAP_FOR_EACH_SAFE (mre, next, struct rank_entry, node, &kw[i].map)
 250         {
 251           hmap_delete (&kw[i].map, &mre->node);
 252           free (mre);
 253         }
 254       hmap_destroy (&kw[i].map);
 255     }
 256
 257   free (kw);
 258 }
 259 \f
 260 static void
 261 show_ranks_box (const struct n_sample_test *nst, const struct kw *kw)
 262 {
 263   struct pivot_table *table = pivot_table_create (N_("Ranks"));
 264
 265   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
 266                           N_("N"), PIVOT_RC_INTEGER,
 267                           N_("Mean Rank"), PIVOT_RC_OTHER);
 268
 269   struct pivot_dimension *variables = pivot_dimension_create (
 270     table, PIVOT_AXIS_ROW, N_("Variables"));
 271
 272   for (size_t i = 0 ; i < nst->n_vars ; ++i)
 273     {
 274       /* Sort the rank entries, by iteratin the hash and putting the entries
 275          into a binary tree. */
 276       struct bt bt = BT_INITIALIZER(compare_rank_entries_3way, nst->vars[i]);
 277       struct rank_entry *re_x;
 278       HMAP_FOR_EACH (re_x, struct rank_entry, node, &kw[i].map)
 279         bt_insert (&bt, &re_x->btn);
 280
 281       /* Report the rank entries in sorted order. */
 282       struct pivot_category *group = pivot_category_create_group__ (
 283         variables->root, pivot_value_new_variable (nst->vars[i]));
 284       int tot = 0;
 285       const struct rank_entry *re;
 286       BT_FOR_EACH (re, struct rank_entry, btn, &bt)
 287         {
 288           struct string str = DS_EMPTY_INITIALIZER;
 289           var_append_value_name (nst->indep_var, &re->group, &str);
 290           int row = pivot_category_create_leaf (
 291             group, pivot_value_new_user_text_nocopy (ds_steal_cstr (&str)));
 292
 293           double entries[] = { re->n, re->sum_of_ranks / re->n };
 294           for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 295             pivot_table_put2 (table, j, row,
 296                               pivot_value_new_number (entries[j]));
 297
 298           tot += re->n;
 299         }
 300
 301       int row = pivot_category_create_leaves (group, N_("Total"));
 302       pivot_table_put2 (table, 0, row, pivot_value_new_number (tot));
 303     }
 304
 305   pivot_table_submit (table);
 306 }
 307
 308 static void
 309 show_sig_box (const struct n_sample_test *nst, const struct kw *kw)
 310 {
 311   struct pivot_table *table = pivot_table_create (N_("Test Statistics"));
 312
 313   pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Statistics"),
 314                           N_("Chi-Square"), PIVOT_RC_OTHER,
 315                           N_("df"), PIVOT_RC_INTEGER,
 316                           N_("Asymp. Sig."), PIVOT_RC_SIGNIFICANCE);
 317
 318   struct pivot_dimension *variables = pivot_dimension_create (
 319     table, PIVOT_AXIS_COLUMN, N_("Variables"));
 320
 321   for (size_t i = 0 ; i < nst->n_vars; ++i)
 322     {
 323       int col = pivot_category_create_leaf (
 324         variables->root, pivot_value_new_variable (nst->vars[i]));
 325
 326       double df = hmap_count (&kw[i].map) - 1;
 327       double sig = gsl_cdf_chisq_Q (kw[i].h, df);
 328       double entries[] = { kw[i].h, df, sig };
 329       for (size_t j = 0; j < sizeof entries / sizeof *entries; j++)
 330         pivot_table_put2 (table, j, col, pivot_value_new_number (entries[j]));
 331     }
 332
 333   pivot_table_submit (table);
 334 }