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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2006, 2007, 2009, 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 #include <config.h>
  18
  19 #include "language/stats/chisquare.h"
  20
  21 #include <gsl/gsl_cdf.h>
  22 #include <math.h>
  23 #include <stdlib.h>
  24
  25 #include "data/case.h"
  26 #include "data/casereader.h"
  27 #include "data/dataset.h"
  28 #include "data/dictionary.h"
  29 #include "data/format.h"
  30 #include "data/value-labels.h"
  31 #include "data/variable.h"
  32 #include "language/stats/freq.h"
  33 #include "language/stats/npar.h"
  34 #include "libpspp/array.h"
  35 #include "libpspp/assertion.h"
  36 #include "libpspp/cast.h"
  37 #include "libpspp/compiler.h"
  38 #include "libpspp/hash-functions.h"
  39 #include "libpspp/message.h"
  40 #include "libpspp/taint.h"
  41 #include "output/tab.h"
  42
  43 #include "gl/xalloc.h"
  44
  45 #include "gettext.h"
  46 #define _(msgid) gettext (msgid)
  47
  48 /* Adds frequency counts of each value of VAR in INPUT between LO and HI to
  49    FREQ_HASH.  LO and HI and each input value is truncated to an integer.
  50    Returns true if successful, false on input error.  It is the caller's
  51    responsibility to initialize FREQ_HASH and to free it when no longer
  52    required, even on failure. */
  53 static bool
  54 create_freq_hash_with_range (const struct dictionary *dict,
  55                              struct casereader *input,
  56                              const struct variable *var,
  57                              double lo_, double hi_,
  58                              struct hmap *freq_hash)
  59 {
  60   struct freq **entries;
  61   bool warn = true;
  62   struct ccase *c;
  63   double lo, hi;
  64   double i_d;
  65
  66   assert (var_is_numeric (var));
  67   lo = trunc (lo_);
  68   hi = trunc (hi_);
  69
  70   /* Populate the hash with zero entries */
  71   entries = xnmalloc (hi - lo + 1, sizeof *entries);
  72   for (i_d = lo; i_d <= hi; i_d += 1.0 )
  73     {
  74       size_t ofs = i_d - lo;
  75       union value value = { i_d };
  76       entries[ofs] = freq_hmap_insert (freq_hash, &value, 0,
  77                                        value_hash (&value, 0, 0));
  78     }
  79
  80   for (; (c = casereader_read (input)) != NULL; case_unref (c))
  81     {
  82       double x = trunc (case_num (c, var));
  83       if (x >= lo && x <= hi)
  84         {
  85           size_t ofs = x - lo;
  86           struct freq *fr = entries[ofs];
  87           fr->count += dict_get_case_weight (dict, c, &warn);
  88         }
  89     }
  90
  91   free (entries);
  92
  93   return casereader_destroy (input);
  94 }
  95
  96 /* Adds frequency counts of each value of VAR in INPUT to FREQ_HASH.  LO and HI
  97    and each input value is truncated to an integer.  Returns true if
  98    successful, false on input error.  It is the caller's responsibility to
  99    initialize FREQ_HASH and to free it when no longer required, even on
 100    failure. */
 101 static bool
 102 create_freq_hash (const struct dictionary *dict,
 103                   struct casereader *input,
 104                   const struct variable *var,
 105                   struct hmap *freq_hash)
 106 {
 107   int width = var_get_width (var);
 108   bool warn = true;
 109   struct ccase *c;
 110
 111   for (; (c = casereader_read (input)) != NULL; case_unref (c))
 112     {
 113       const union value *value = case_data (c, var);
 114       size_t hash = value_hash (value, width, 0);
 115       double weight = dict_get_case_weight (dict, c, &warn);
 116       struct freq *f;
 117
 118       f = freq_hmap_search (freq_hash, value, width, hash);
 119       if (f == NULL)
 120         f = freq_hmap_insert (freq_hash, value, width, hash);
 121
 122       f->count += weight;
 123     }
 124
 125   return casereader_destroy (input);
 126 }
 127
 128 static struct tab_table *
 129 create_variable_frequency_table (const struct dictionary *dict,
 130                                  struct casereader *input,
 131                                  const struct chisquare_test *test,
 132                                  int v, struct hmap *freq_hash)
 133
 134 {
 135   int i;
 136   const struct one_sample_test *ost = (const struct one_sample_test*)test;
 137   int n_cells;
 138   struct tab_table *table ;
 139   const struct variable *var =  ost->vars[v];
 140
 141   const struct variable *wvar = dict_get_weight (dict);
 142   const struct fmt_spec *wfmt = wvar ? var_get_print_format (wvar) : & F_8_0;
 143
 144   hmap_init (freq_hash);
 145   if (!create_freq_hash (dict, input, var, freq_hash))
 146     {
 147       freq_hmap_destroy (freq_hash, var_get_width (var));
 148       return NULL;
 149     }
 150
 151   n_cells = hmap_count (freq_hash);
 152
 153   if ( test->n_expected > 0 && n_cells != test->n_expected )
 154     {
 155       msg(ME, _("CHISQUARE test specified %d expected values, but"
 156                 " %d distinct values were encountered in variable %s."),
 157           test->n_expected, n_cells,
 158           var_get_name (var)
 159           );
 160       freq_hmap_destroy (freq_hash, var_get_width (var));
 161       return NULL;
 162     }
 163
 164   table = tab_create(4, n_cells + 2);
 165   tab_set_format (table, RC_WEIGHT, wfmt);
 166
 167   tab_title (table, "%s", var_to_string(var));
 168   tab_text (table, 1, 0, TAB_LEFT, _("Observed N"));
 169   tab_text (table, 2, 0, TAB_LEFT, _("Expected N"));
 170   tab_text (table, 3, 0, TAB_LEFT, _("Residual"));
 171
 172   tab_headers (table, 1, 0, 1, 0);
 173
 174   tab_box (table, TAL_1, TAL_1, -1, -1,
 175            0, 0, tab_nc (table) - 1, tab_nr(table) - 1 );
 176
 177   tab_hline (table, TAL_1, 0, tab_nc(table) - 1, 1);
 178
 179   tab_vline (table, TAL_2, 1, 0, tab_nr(table) - 1);
 180   for ( i = 2 ; i < 4 ; ++i )
 181     tab_vline (table, TAL_1, i, 0, tab_nr(table) - 1);
 182
 183
 184   tab_text (table, 0, tab_nr (table) - 1, TAB_LEFT, _("Total"));
 185
 186   return table;
 187 }
 188
 189
 190 static struct tab_table *
 191 create_combo_frequency_table (const struct dictionary *dict, const struct chisquare_test *test)
 192 {
 193   int i;
 194   const struct one_sample_test *ost = (const struct one_sample_test*)test;
 195
 196   struct tab_table *table ;
 197
 198   const struct variable *wvar = dict_get_weight (dict);
 199   const struct fmt_spec *wfmt = wvar ? var_get_print_format (wvar) : & F_8_0;
 200
 201   int n_cells = test->hi - test->lo + 1;
 202
 203   table = tab_create(1 + ost->n_vars * 4, n_cells + 3);
 204   tab_set_format (table, RC_WEIGHT, wfmt);
 205
 206   tab_title (table, _("Frequencies"));
 207   for ( i = 0 ; i < ost->n_vars ; ++i )
 208     {
 209       const struct variable *var = ost->vars[i];
 210       tab_text (table, i * 4 + 1, 1, TAB_LEFT, _("Category"));
 211       tab_text (table, i * 4 + 2, 1, TAB_LEFT, _("Observed N"));
 212       tab_text (table, i * 4 + 3, 1, TAB_LEFT, _("Expected N"));
 213       tab_text (table, i * 4 + 4, 1, TAB_LEFT, _("Residual"));
 214
 215       tab_vline (table, TAL_2, i * 4 + 1,
 216                  0, tab_nr (table) - 1);
 217
 218       tab_vline (table, TAL_1, i * 4 + 2,
 219                  0, tab_nr (table) - 1);
 220
 221       tab_vline (table, TAL_1, i * 4 + 3,
 222                  1, tab_nr (table) - 1);
 223
 224       tab_vline (table, TAL_1, i * 4 + 4,
 225                  1, tab_nr (table) - 1);
 226
 227
 228       tab_joint_text (table,
 229                       i * 4 + 1, 0,
 230                       i * 4 + 4, 0,
 231                       TAB_CENTER,
 232                       var_to_string (var));
 233     }
 234
 235   for ( i = test->lo ; i <= test->hi ; ++i )
 236     tab_double (table, 0, 2 + i - test->lo, TAB_LEFT, 1 + i - test->lo, NULL, RC_INTEGER);
 237
 238   tab_headers (table, 1, 0, 2, 0);
 239
 240   tab_box (table, TAL_1, TAL_1, -1, -1,
 241            0, 0, tab_nc (table) - 1, tab_nr(table) - 1 );
 242
 243   tab_hline (table, TAL_1, 1, tab_nc(table) - 1, 1);
 244   tab_hline (table, TAL_1, 0, tab_nc(table) - 1, 2);
 245
 246   tab_text (table, 0, tab_nr (table) - 1, TAB_LEFT, _("Total"));
 247
 248   return table;
 249 }
 250
 251
 252 static struct tab_table *
 253 create_stats_table (const struct chisquare_test *test)
 254 {
 255   const struct one_sample_test *ost = (const struct one_sample_test*) test;
 256
 257   struct tab_table *table;
 258   table = tab_create (1 + ost->n_vars, 4);
 259   tab_title (table, _("Test Statistics"));
 260   tab_headers (table, 1, 0, 1, 0);
 261
 262   tab_box (table, TAL_1, TAL_1, -1, -1,
 263            0, 0, tab_nc(table) - 1, tab_nr(table) - 1 );
 264
 265   tab_box (table, -1, -1, -1, TAL_1,
 266            1, 0, tab_nc(table) - 1, tab_nr(table) - 1 );
 267
 268
 269   tab_vline (table, TAL_2, 1, 0, tab_nr (table) - 1);
 270   tab_hline (table, TAL_1, 0, tab_nc (table) - 1, 1);
 271
 272
 273   tab_text (table, 0, 1, TAB_LEFT, _("Chi-Square"));
 274   tab_text (table, 0, 2, TAB_LEFT, _("df"));
 275   tab_text (table, 0, 3, TAB_LEFT, _("Asymp. Sig."));
 276
 277   return table;
 278 }
 279
 280
 281 void
 282 chisquare_execute (const struct dataset *ds,
 283                    struct casereader *input,
 284                    enum mv_class exclude,
 285                    const struct npar_test *test,
 286                    bool exact UNUSED,
 287                    double timer UNUSED)
 288 {
 289   const struct dictionary *dict = dataset_dict (ds);
 290   int v, i;
 291   struct chisquare_test *cst = UP_CAST (test, struct chisquare_test,
 292                                         parent.parent);
 293   struct one_sample_test *ost = &cst->parent;
 294   int n_cells = 0;
 295   double total_expected = 0.0;
 296
 297   double *df = xzalloc (sizeof (*df) * ost->n_vars);
 298   double *xsq = xzalloc (sizeof (*df) * ost->n_vars);
 299   bool ok;
 300
 301   for ( i = 0 ; i < cst->n_expected ; ++i )
 302     total_expected += cst->expected[i];
 303
 304   if ( cst->ranged == false )
 305     {
 306       for ( v = 0 ; v < ost->n_vars ; ++v )
 307         {
 308           const struct variable *var = ost->vars[v];
 309           double total_obs = 0.0;
 310           struct hmap freq_hash;
 311           struct casereader *reader =
 312             casereader_create_filter_missing (casereader_clone (input),
 313                                               &var, 1, exclude,
 314                                               NULL, NULL);
 315           struct tab_table *freq_table =
 316             create_variable_frequency_table (dict, reader, cst, v, &freq_hash);
 317
 318           struct freq **ff;
 319
 320           if ( NULL == freq_table )
 321             continue;
 322           ff = freq_hmap_sort (&freq_hash, var_get_width (var));
 323
 324           n_cells = hmap_count (&freq_hash);
 325
 326           for ( i = 0 ; i < n_cells ; ++i )
 327             total_obs += ff[i]->count;
 328
 329           xsq[v] = 0.0;
 330           for ( i = 0 ; i < n_cells ; ++i )
 331             {
 332               struct string str;
 333               double exp;
 334               const union value *observed_value = &ff[i]->values[0];
 335
 336               ds_init_empty (&str);
 337               var_append_value_name (var, observed_value, &str);
 338
 339               /* The key */
 340               tab_text (freq_table, 0, i + 1, TAB_LEFT, ds_cstr (&str));
 341               ds_destroy (&str);
 342
 343
 344               /* The observed N */
 345               tab_double (freq_table, 1, i + 1, TAB_NONE,
 346                           ff[i]->count, NULL, RC_WEIGHT);
 347
 348               if ( cst->n_expected > 0 )
 349                 exp = cst->expected[i] * total_obs / total_expected ;
 350               else
 351                 exp = total_obs / (double) n_cells;
 352
 353               tab_double (freq_table, 2, i + 1, TAB_NONE,
 354                           exp, NULL, RC_OTHER);
 355
 356               /* The residual */
 357               tab_double (freq_table, 3, i + 1, TAB_NONE,
 358                           ff[i]->count - exp, NULL, RC_OTHER);
 359
 360               xsq[v] += (ff[i]->count - exp) * (ff[i]->count - exp) / exp;
 361             }
 362
 363           df[v] = n_cells - 1.0;
 364
 365           tab_double (freq_table, 1, i + 1, TAB_NONE,
 366                       total_obs, NULL, RC_WEIGHT);
 367
 368           tab_submit (freq_table);
 369
 370           freq_hmap_destroy (&freq_hash, var_get_width (var));
 371           free (ff);
 372         }
 373     }
 374   else  /* ranged == true */
 375     {
 376       struct tab_table *freq_table = create_combo_frequency_table (dict, cst);
 377
 378       n_cells = cst->hi - cst->lo + 1;
 379
 380       for ( v = 0 ; v < ost->n_vars ; ++v )
 381         {
 382           const struct variable *var = ost->vars[v];
 383           double total_obs = 0.0;
 384           struct casereader *reader =
 385             casereader_create_filter_missing (casereader_clone (input),
 386                                               &var, 1, exclude,
 387                                               NULL, NULL);
 388           struct hmap freq_hash;
 389           struct freq **ff;
 390
 391           hmap_init (&freq_hash);
 392           if (!create_freq_hash_with_range (dict, reader, var,
 393                                             cst->lo, cst->hi, &freq_hash))
 394             {
 395               freq_hmap_destroy (&freq_hash, var_get_width (var));
 396               continue;
 397             }
 398
 399           ff = freq_hmap_sort (&freq_hash, var_get_width (var));
 400
 401           for ( i = 0 ; i < hmap_count (&freq_hash) ; ++i )
 402             total_obs += ff[i]->count;
 403
 404           xsq[v] = 0.0;
 405           for ( i = 0 ; i < hmap_count (&freq_hash) ; ++i )
 406             {
 407               struct string str;
 408               double exp;
 409
 410               const union value *observed_value = &ff[i]->values[0];
 411
 412               ds_init_empty (&str);
 413               var_append_value_name (ost->vars[v], observed_value, &str);
 414               /* The key */
 415               tab_text  (freq_table, v * 4 + 1, i + 2 , TAB_LEFT,
 416                          ds_cstr (&str));
 417               ds_destroy (&str);
 418
 419               /* The observed N */
 420               tab_double (freq_table, v * 4 + 2, i + 2 , TAB_NONE,
 421                           ff[i]->count, NULL, RC_WEIGHT);
 422
 423               if ( cst->n_expected > 0 )
 424                 exp = cst->expected[i] * total_obs / total_expected ;
 425               else
 426                 exp = total_obs / (double) hmap_count (&freq_hash);
 427
 428               /* The expected N */
 429               tab_double (freq_table, v * 4 + 3, i + 2 , TAB_NONE,
 430                           exp, NULL, RC_OTHER);
 431
 432               /* The residual */
 433               tab_double (freq_table, v * 4 + 4, i + 2 , TAB_NONE,
 434                           ff[i]->count - exp, NULL, RC_OTHER);
 435
 436               xsq[v] += (ff[i]->count - exp) * (ff[i]->count - exp) / exp;
 437             }
 438
 439
 440           tab_double (freq_table, v * 4 + 2, tab_nr (freq_table) - 1, TAB_NONE,
 441                       total_obs, NULL, RC_WEIGHT);
 442
 443           df[v] = n_cells - 1.0;
 444
 445           freq_hmap_destroy (&freq_hash, var_get_width (var));
 446           free (ff);
 447         }
 448
 449       tab_submit (freq_table);
 450     }
 451   ok = !taint_has_tainted_successor (casereader_get_taint (input));
 452   casereader_destroy (input);
 453
 454   if (ok)
 455     {
 456       struct tab_table *stats_table = create_stats_table (cst);
 457
 458       /* Populate the summary statistics table */
 459       for ( v = 0 ; v < ost->n_vars ; ++v )
 460         {
 461           const struct variable *var = ost->vars[v];
 462
 463           tab_text (stats_table, 1 + v, 0, TAB_CENTER, var_get_name (var));
 464
 465           tab_double (stats_table, 1 + v, 1, TAB_NONE, xsq[v], NULL, RC_OTHER);
 466           tab_double (stats_table, 1 + v, 2, TAB_NONE, df[v], NULL, RC_INTEGER);
 467
 468           tab_double (stats_table, 1 + v, 3, TAB_NONE,
 469                       gsl_cdf_chisq_Q (xsq[v], df[v]), NULL, RC_PVALUE);
 470         }
 471       tab_submit (stats_table);
 472     }
 473
 474   free (xsq);
 475   free (df);
 476 }
 477