src/language/stats/chisquare.c

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