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