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               continue;
 201             }
 202
 203           if (n_sample_test)
 204             {
 205               int width = var_get_width (nst->indep_var);
 206               /* Ignore out of range values */
 207               if (
 208                   value_compare_3way (indep_val, &nst->val1, width) < 0
 209                 ||
 210                   value_compare_3way (indep_val, &nst->val2, width) > 0
 211                    )
 212                 {
 213                   continue;
 214                 }
 215             }
 216
 217           vn = find_value (&map, indep_val, nst->indep_var);
 218           if ( vn == NULL)
 219             {
 220               if ( n_sample_test == true)
 221                 {
 222                   int width = var_get_width (nst->indep_var);
 223                   vn = xzalloc (sizeof *vn);
 224                   value_clone (&vn->val,  indep_val, width);
 225
 226                   hmap_insert (&map, &vn->node, value_hash (indep_val, width, 0));
 227                 }
 228               else
 229                 {
 230                   continue;
 231                 }
 232             }
 233
 234           if (val->f <= median)
 235             vn->le += weight;
 236           else
 237             vn->gt += weight;
 238
 239           count += weight;
 240         }
 241       casereader_destroy (r);
 242
 243       {
 244         int x = 0;
 245         struct val_node *vn = NULL;
 246         double r_0 = 0;
 247         double r_1 = 0;
 248         HMAP_FOR_EACH (vn, struct val_node, node, &map)
 249           {
 250             r_0 += vn->le;
 251             r_1 += vn->gt;
 252           }
 253
 254         results[v].n = count;
 255         results[v].sorted_array = xcalloc (hmap_count (&map), sizeof (void*));
 256         results[v].var = var;
 257
 258         HMAP_FOR_EACH (vn, struct val_node, node, &map)
 259           {
 260             double e_0j = r_0 * (vn->le + vn->gt) / count;
 261             double e_1j = r_1 * (vn->le + vn->gt) / count;
 262
 263             results[v].chisq += pow2 (vn->le - e_0j) / e_0j;
 264             results[v].chisq += pow2 (vn->gt - e_1j) / e_1j;
 265
 266             results[v].sorted_array[x++] = vn;
 267           }
 268
 269         n_vals = x;
 270         hmap_destroy (&map);
 271
 272         sort (results[v].sorted_array, x, sizeof (*results[v].sorted_array),
 273               val_node_cmp_3way, nst->indep_var);
 274
 275       }
 276     }
 277
 278   casereader_destroy (input);
 279
 280   show_frequencies (nst, results,  n_vals, dict);
 281   show_test_statistics (nst, results, n_vals, dict);
 282
 283   for (v = 0; v < nst->n_vars; ++v)
 284     {
 285       int i;
 286       const struct results *rs = results + v;
 287
 288       for (i = 0; i < n_vals; ++i)
 289         {
 290           struct val_node *vn = rs->sorted_array[i];
 291           value_destroy (&vn->val, var_get_width (nst->indep_var));
 292           free (vn);
 293         }
 294       free (rs->sorted_array);
 295     }
 296   free (results);
 297 }
 298
 299
 300
 301 static void
 302 show_frequencies (const struct n_sample_test *nst, const struct results *results,  int n_vals, const struct dictionary *dict)
 303 {
 304   const struct variable *weight = dict_get_weight (dict);
 305   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 306
 307   int i;
 308   int v;
 309
 310   const int row_headers = 2;
 311   const int column_headers = 2;
 312   const int nc = row_headers + n_vals;
 313   const int nr = column_headers + nst->n_vars * 2;
 314
 315   struct tab_table *table = tab_create (nc, nr);
 316
 317   tab_headers (table, row_headers, 0, column_headers, 0);
 318
 319   tab_title (table, _("Frequencies"));
 320
 321   /* Box around the table and vertical lines inside*/
 322   tab_box (table, TAL_2, TAL_2, -1, TAL_1,
 323            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 324
 325   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 326   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 327
 328   tab_joint_text (table,
 329                   row_headers, 0, row_headers + n_vals - 1, 0,
 330                   TAT_TITLE | TAB_CENTER, var_to_string (nst->indep_var));
 331
 332
 333   tab_hline (table, TAL_1, row_headers, tab_nc (table) - 1, 1);
 334
 335
 336   for (i = 0; i < n_vals; ++i)
 337     {
 338       const struct results *rs = results + 0;
 339       struct string label;
 340       ds_init_empty (&label);
 341
 342       var_append_value_name (nst->indep_var, &rs->sorted_array[i]->val,
 343                             &label);
 344
 345       tab_text (table, row_headers + i, 1,
 346                 TAT_TITLE | TAB_LEFT, ds_cstr (&label));
 347
 348       ds_destroy (&label);
 349     }
 350
 351   for (v = 0; v < nst->n_vars; ++v)
 352     {
 353       const struct results *rs = &results[v];
 354       tab_text (table,  0, column_headers + v * 2,
 355                 TAT_TITLE | TAB_LEFT, var_to_string (rs->var) );
 356
 357       tab_text (table,  1, column_headers + v * 2,
 358                 TAT_TITLE | TAB_LEFT, _("> Median") );
 359
 360       tab_text (table,  1, column_headers + v * 2 + 1,
 361                 TAT_TITLE | TAB_LEFT, _("≤ Median") );
 362
 363       if ( v > 0)
 364         tab_hline (table, TAL_1, 0, tab_nc (table) - 1, column_headers + v * 2);
 365     }
 366
 367   for (v = 0; v < nst->n_vars; ++v)
 368     {
 369       int i;
 370       const struct results *rs = &results[v];
 371
 372       for (i = 0; i < n_vals; ++i)
 373         {
 374           const struct val_node *vn = rs->sorted_array[i];
 375           tab_double (table, row_headers + i, column_headers + v * 2,
 376                       0, vn->gt, wfmt);
 377
 378           tab_double (table, row_headers + i, column_headers + v * 2 + 1,
 379                     0, vn->le, wfmt);
 380         }
 381     }
 382
 383   tab_submit (table);
 384 }
 385
 386
 387 static void
 388 show_test_statistics (const struct n_sample_test *nst,
 389                       const struct results *results,
 390                       int n_vals,
 391                       const struct dictionary *dict)
 392 {
 393   const struct variable *weight = dict_get_weight (dict);
 394   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 395
 396   int v;
 397
 398   const int row_headers = 1;
 399   const int column_headers = 1;
 400   const int nc = row_headers + 5;
 401   const int nr = column_headers + nst->n_vars;
 402
 403   struct tab_table *table = tab_create (nc, nr);
 404
 405   tab_headers (table, row_headers, 0, column_headers, 0);
 406
 407   tab_title (table, _("Test Statistics"));
 408
 409
 410   tab_box (table, TAL_2, TAL_2, -1, TAL_1,
 411            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 412
 413   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 414   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 415
 416   tab_text (table, row_headers + 0, 0,
 417             TAT_TITLE | TAB_CENTER, _("N"));
 418
 419   tab_text (table, row_headers + 1, 0,
 420             TAT_TITLE | TAB_CENTER, _("Median"));
 421
 422   tab_text (table, row_headers + 2, 0,
 423             TAT_TITLE | TAB_CENTER, _("Chi-Square"));
 424
 425   tab_text (table, row_headers + 3, 0,
 426             TAT_TITLE | TAB_CENTER, _("df"));
 427
 428   tab_text (table, row_headers + 4, 0,
 429             TAT_TITLE | TAB_CENTER, _("Asymp. Sig."));
 430
 431
 432   for (v = 0; v < nst->n_vars; ++v)
 433     {
 434       double df = n_vals - 1;
 435       const struct results *rs = &results[v];
 436       tab_text (table,  0, column_headers + v,
 437                 TAT_TITLE | TAB_LEFT, var_to_string (rs->var));
 438
 439
 440       tab_double (table, row_headers + 0, column_headers + v,
 441                   0, rs->n, wfmt);
 442
 443       tab_double (table, row_headers + 1, column_headers + v,
 444                   0, rs->median, NULL);
 445
 446       tab_double (table, row_headers + 2, column_headers + v,
 447                   0, rs->chisq, NULL);
 448
 449       tab_double (table, row_headers + 3, column_headers + v,
 450                   0, df, wfmt);
 451
 452       tab_double (table, row_headers + 4, column_headers + v,
 453                   0, gsl_cdf_chisq_Q (rs->chisq, df), NULL);
 454     }
 455
 456   tab_submit (table);
 457 }