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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 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 #include "median.h"
  20
  21 #include <gsl/gsl_cdf.h>
  22
  23 #include "data/format.h"
  24
  25
  26 #include "data/variable.h"
  27 #include "data/case.h"
  28 #include "data/dictionary.h"
  29 #include "data/dataset.h"
  30 #include "data/casereader.h"
  31 #include "data/casewriter.h"
  32 #include "data/subcase.h"
  33 #include "data/value.h"
  34
  35 #include "math/percentiles.h"
  36 #include "math/sort.h"
  37
  38 #include "libpspp/cast.h"
  39 #include "libpspp/hmap.h"
  40 #include "libpspp/array.h"
  41 #include "libpspp/str.h"
  42 #include "libpspp/misc.h"
  43
  44 #include "output/tab.h"
  45
  46 #include "gettext.h"
  47 #define _(msgid) gettext (msgid)
  48
  49
  50 struct val_node
  51 {
  52   struct hmap_node node;
  53   union value val;
  54   casenumber le;
  55   casenumber gt;
  56 };
  57
  58 struct results
  59 {
  60   const struct variable *var;
  61   struct val_node **sorted_array;
  62   double n;
  63   double median;
  64   double chisq;
  65 };
  66
  67
  68
  69 static int
  70 val_node_cmp_3way (const void *a_, const void *b_, const void *aux)
  71 {
  72   const struct variable *indep_var = aux;
  73   const struct val_node *const *a = a_;
  74   const struct val_node *const *b = b_;
  75
  76   return value_compare_3way (&(*a)->val, &(*b)->val, var_get_width (indep_var));
  77 }
  78
  79 static void
  80 show_frequencies (const struct n_sample_test *nst, const struct results *results,  int n_vals, const struct dictionary *);
  81
  82 static void
  83 show_test_statistics (const struct n_sample_test *nst, const struct results *results, int, const struct dictionary *);
  84
  85
  86 static struct val_node *
  87 find_value (const struct hmap *map, const union value *val,
  88             const struct variable *var)
  89 {
  90   struct val_node *foo = NULL;
  91   size_t hash = value_hash (val, var_get_width (var), 0);
  92   HMAP_FOR_EACH_WITH_HASH (foo, struct val_node, node, hash, map)
  93     if (value_equal (val, &foo->val, var_get_width (var)))
  94       break;
  95
  96   return foo;
  97 }
  98
  99 void
 100 median_execute (const struct dataset *ds,
 101                 struct casereader *input,
 102                 enum mv_class exclude,
 103                 const struct npar_test *test,
 104                 bool exact UNUSED,
 105                 double timer UNUSED)
 106 {
 107   const struct dictionary *dict = dataset_dict (ds);
 108   const struct variable *wvar = dict_get_weight (dict);
 109   bool warn = true;
 110   int v;
 111   const struct median_test *mt = UP_CAST (test, const struct median_test,
 112                                           parent.parent);
 113
 114   const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test,
 115                                           parent);
 116
 117   const bool n_sample_test = (value_compare_3way (&nst->val2, &nst->val1,
 118                                        var_get_width (nst->indep_var)) > 0);
 119
 120   struct results *results = xcalloc (nst->n_vars, sizeof (*results));
 121   int n_vals = 0;
 122   for (v = 0; v < nst->n_vars; ++v)
 123     {
 124       double count = 0;
 125       double cc = 0;
 126       double median = mt->median;
 127       const struct variable *var = nst->vars[v];
 128       struct ccase *c;
 129       struct hmap map = HMAP_INITIALIZER (map);
 130       struct casereader *r = casereader_clone (input);
 131
 132
 133
 134       if (n_sample_test == false)
 135         {
 136           struct val_node *vn = xzalloc (sizeof *vn);
 137           value_clone (&vn->val,  &nst->val1, var_get_width (nst->indep_var));
 138           hmap_insert (&map, &vn->node, value_hash (&nst->val1,
 139                                             var_get_width (nst->indep_var), 0));
 140
 141           vn = xzalloc (sizeof *vn);
 142           value_clone (&vn->val,  &nst->val2, var_get_width (nst->indep_var));
 143           hmap_insert (&map, &vn->node, value_hash (&nst->val2,
 144                                             var_get_width (nst->indep_var), 0));
 145         }
 146
 147       if ( median == SYSMIS)
 148         {
 149           struct percentile *ptl;
 150           struct order_stats *os;
 151
 152           struct casereader *rr;
 153           struct subcase sc;
 154           struct casewriter *writer;
 155           subcase_init_var (&sc, var, SC_ASCEND);
 156           rr = casereader_clone (r);
 157           writer = sort_create_writer (&sc, casereader_get_proto (rr));
 158
 159           for (; (c = casereader_read (rr)) != NULL; )
 160             {
 161               if ( var_is_value_missing (var, case_data (c, var), exclude))
 162                 {
 163                   case_unref (c);
 164                   continue;
 165                 }
 166
 167               cc += dict_get_case_weight (dict, c, &warn);
 168               casewriter_write (writer, c);
 169             }
 170           subcase_destroy (&sc);
 171           casereader_destroy (rr);
 172
 173           rr = casewriter_make_reader (writer);
 174
 175           ptl = percentile_create (0.5, cc);
 176           os = &ptl->parent;
 177
 178           order_stats_accumulate (&os, 1,
 179                                   rr,
 180                                   wvar,
 181                                   var,
 182                                   exclude);
 183
 184           median = percentile_calculate (ptl, PC_HAVERAGE);
 185           statistic_destroy (&ptl->parent.parent);
 186         }
 187
 188       results[v].median = median;
 189
 190
 191       for (; (c = casereader_read (r)) != NULL; case_unref (c))
 192         {
 193           struct val_node *vn ;
 194           const double weight = dict_get_case_weight (dict, c, &warn);
 195           const union value *val = case_data (c, var);
 196           const union value *indep_val = case_data (c, nst->indep_var);
 197
 198           if ( var_is_value_missing (var, case_data (c, var), exclude))
 199             {
 200               case_unref (c);
 201               continue;
 202             }
 203
 204           if (n_sample_test)
 205             {
 206               int width = var_get_width (nst->indep_var);
 207               /* Ignore out of range values */
 208               if (
 209                   value_compare_3way (indep_val, &nst->val1, width) < 0
 210                 ||
 211                   value_compare_3way (indep_val, &nst->val2, width) > 0
 212                    )
 213                 {
 214                   case_unref (c);
 215                   continue;
 216                 }
 217             }
 218
 219           vn = find_value (&map, indep_val, nst->indep_var);
 220           if ( vn == NULL)
 221             {
 222               if ( n_sample_test == true)
 223                 {
 224                   int width = var_get_width (nst->indep_var);
 225                   vn = xzalloc (sizeof *vn);
 226                   value_clone (&vn->val,  indep_val, width);
 227
 228                   hmap_insert (&map, &vn->node, value_hash (indep_val, width, 0));
 229                 }
 230               else
 231                 {
 232                   continue;
 233                 }
 234             }
 235
 236           if (val->f <= median)
 237             vn->le += weight;
 238           else
 239             vn->gt += weight;
 240
 241           count += weight;
 242         }
 243       casereader_destroy (r);
 244
 245       {
 246         int x = 0;
 247         struct val_node *vn = NULL;
 248         double r_0 = 0;
 249         double r_1 = 0;
 250         HMAP_FOR_EACH (vn, struct val_node, node, &map)
 251           {
 252             r_0 += vn->le;
 253             r_1 += vn->gt;
 254           }
 255
 256         results[v].n = count;
 257         results[v].sorted_array = xcalloc (hmap_count (&map), sizeof (void*));
 258         results[v].var = var;
 259
 260         HMAP_FOR_EACH (vn, struct val_node, node, &map)
 261           {
 262             double e_0j = r_0 * (vn->le + vn->gt) / count;
 263             double e_1j = r_1 * (vn->le + vn->gt) / count;
 264
 265             results[v].chisq += pow2 (vn->le - e_0j) / e_0j;
 266             results[v].chisq += pow2 (vn->gt - e_1j) / e_1j;
 267
 268             results[v].sorted_array[x++] = vn;
 269           }
 270
 271         n_vals = x;
 272         hmap_destroy (&map);
 273
 274         sort (results[v].sorted_array, x, sizeof (*results[v].sorted_array),
 275               val_node_cmp_3way, nst->indep_var);
 276
 277       }
 278     }
 279
 280   casereader_destroy (input);
 281
 282   show_frequencies (nst, results,  n_vals, dict);
 283   show_test_statistics (nst, results, n_vals, dict);
 284
 285   for (v = 0; v < nst->n_vars; ++v)
 286     {
 287       int i;
 288       const struct results *rs = results + v;
 289
 290       for (i = 0; i < n_vals; ++i)
 291         {
 292           struct val_node *vn = rs->sorted_array[i];
 293           value_destroy (&vn->val, var_get_width (nst->indep_var));
 294           free (vn);
 295         }
 296       free (rs->sorted_array);
 297     }
 298   free (results);
 299 }
 300
 301
 302
 303 static void
 304 show_frequencies (const struct n_sample_test *nst, const struct results *results,  int n_vals, const struct dictionary *dict)
 305 {
 306   const struct variable *weight = dict_get_weight (dict);
 307   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 308
 309   int i;
 310   int v;
 311
 312   const int row_headers = 2;
 313   const int column_headers = 2;
 314   const int nc = row_headers + n_vals;
 315   const int nr = column_headers + nst->n_vars * 2;
 316
 317   struct tab_table *table = tab_create (nc, nr);
 318
 319   tab_headers (table, row_headers, 0, column_headers, 0);
 320
 321   tab_title (table, _("Frequencies"));
 322
 323   /* Box around the table and vertical lines inside*/
 324   tab_box (table, TAL_2, TAL_2, -1, TAL_1,
 325            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 326
 327   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 328   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 329
 330   tab_joint_text (table,
 331                   row_headers, 0, row_headers + n_vals - 1, 0,
 332                   TAT_TITLE | TAB_CENTER, var_to_string (nst->indep_var));
 333
 334
 335   tab_hline (table, TAL_1, row_headers, tab_nc (table) - 1, 1);
 336
 337
 338   for (i = 0; i < n_vals; ++i)
 339     {
 340       const struct results *rs = results + 0;
 341       struct string label;
 342       ds_init_empty (&label);
 343
 344       var_append_value_name (nst->indep_var, &rs->sorted_array[i]->val,
 345                             &label);
 346
 347       tab_text (table, row_headers + i, 1,
 348                 TAT_TITLE | TAB_LEFT, ds_cstr (&label));
 349
 350       ds_destroy (&label);
 351     }
 352
 353   for (v = 0; v < nst->n_vars; ++v)
 354     {
 355       const struct results *rs = &results[v];
 356       tab_text (table,  0, column_headers + v * 2,
 357                 TAT_TITLE | TAB_LEFT, var_to_string (rs->var) );
 358
 359       tab_text (table,  1, column_headers + v * 2,
 360                 TAT_TITLE | TAB_LEFT, _("> Median") );
 361
 362       tab_text (table,  1, column_headers + v * 2 + 1,
 363                 TAT_TITLE | TAB_LEFT, _("≤ Median") );
 364
 365       if ( v > 0)
 366         tab_hline (table, TAL_1, 0, tab_nc (table) - 1, column_headers + v * 2);
 367     }
 368
 369   for (v = 0; v < nst->n_vars; ++v)
 370     {
 371       int i;
 372       const struct results *rs = &results[v];
 373
 374       for (i = 0; i < n_vals; ++i)
 375         {
 376           const struct val_node *vn = rs->sorted_array[i];
 377           tab_double (table, row_headers + i, column_headers + v * 2,
 378                       0, vn->gt, wfmt);
 379
 380           tab_double (table, row_headers + i, column_headers + v * 2 + 1,
 381                     0, vn->le, wfmt);
 382         }
 383     }
 384
 385   tab_submit (table);
 386 }
 387
 388
 389 static void
 390 show_test_statistics (const struct n_sample_test *nst,
 391                       const struct results *results,
 392                       int n_vals,
 393                       const struct dictionary *dict)
 394 {
 395   const struct variable *weight = dict_get_weight (dict);
 396   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 397
 398   int v;
 399
 400   const int row_headers = 1;
 401   const int column_headers = 1;
 402   const int nc = row_headers + 5;
 403   const int nr = column_headers + nst->n_vars;
 404
 405   struct tab_table *table = tab_create (nc, nr);
 406
 407   tab_headers (table, row_headers, 0, column_headers, 0);
 408
 409   tab_title (table, _("Test Statistics"));
 410
 411
 412   tab_box (table, TAL_2, TAL_2, -1, TAL_1,
 413            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 414
 415   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 416   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 417
 418   tab_text (table, row_headers + 0, 0,
 419             TAT_TITLE | TAB_CENTER, _("N"));
 420
 421   tab_text (table, row_headers + 1, 0,
 422             TAT_TITLE | TAB_CENTER, _("Median"));
 423
 424   tab_text (table, row_headers + 2, 0,
 425             TAT_TITLE | TAB_CENTER, _("Chi-Square"));
 426
 427   tab_text (table, row_headers + 3, 0,
 428             TAT_TITLE | TAB_CENTER, _("df"));
 429
 430   tab_text (table, row_headers + 4, 0,
 431             TAT_TITLE | TAB_CENTER, _("Asymp. Sig."));
 432
 433
 434   for (v = 0; v < nst->n_vars; ++v)
 435     {
 436       double df = n_vals - 1;
 437       const struct results *rs = &results[v];
 438       tab_text (table,  0, column_headers + v,
 439                 TAT_TITLE | TAB_LEFT, var_to_string (rs->var));
 440
 441
 442       tab_double (table, row_headers + 0, column_headers + v,
 443                   0, rs->n, wfmt);
 444
 445       tab_double (table, row_headers + 1, column_headers + v,
 446                   0, rs->median, NULL);
 447
 448       tab_double (table, row_headers + 2, column_headers + v,
 449                   0, rs->chisq, NULL);
 450
 451       tab_double (table, row_headers + 3, column_headers + v,
 452                   0, df, wfmt);
 453
 454       tab_double (table, row_headers + 4, column_headers + v,
 455                   0, gsl_cdf_chisq_Q (rs->chisq, df), NULL);
 456     }
 457
 458   tab_submit (table);
 459 }