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