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