src/language/stats/chisquare.c

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