src/language/stats/chisquare.c

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