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