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