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