Write message to status bar on saving syntax files
[pspp-builds.git] / 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                    bool exact UNUSED,
325                    double timer UNUSED)
326 {
327   const struct dictionary *dict = dataset_dict (ds);
328   int v, i;
329   struct one_sample_test *ost = (struct one_sample_test *) test;
330   struct chisquare_test *cst = (struct chisquare_test *) test;
331   int n_cells = 0;
332   double total_expected = 0.0;
333
334   double *df = xzalloc (sizeof (*df) * ost->n_vars);
335   double *xsq = xzalloc (sizeof (*df) * ost->n_vars);
336   bool ok;
337
338   for ( i = 0 ; i < cst->n_expected ; ++i )
339     total_expected += cst->expected[i];
340
341   if ( cst->ranged == false )
342     {
343       for ( v = 0 ; v < ost->n_vars ; ++v )
344         {
345           double total_obs = 0.0;
346           struct hsh_table *freq_hash = NULL;
347           struct casereader *reader =
348             casereader_create_filter_missing (casereader_clone (input),
349                                               &ost->vars[v], 1, exclude,
350                                               NULL, NULL);
351           struct tab_table *freq_table =
352             create_variable_frequency_table(dict, reader, cst, v, &freq_hash);
353
354           struct freq **ff;
355
356           if ( NULL == freq_table )
357             continue;
358           ff = (struct freq **) hsh_sort (freq_hash);
359
360           n_cells = hsh_count (freq_hash);
361
362           for ( i = 0 ; i < n_cells ; ++i )
363             total_obs += ff[i]->count;
364
365           xsq[v] = 0.0;
366           for ( i = 0 ; i < n_cells ; ++i )
367             {
368               struct string str;
369               double exp;
370               const union value *observed_value = ff[i]->value;
371
372               ds_init_empty (&str);
373               var_append_value_name (ost->vars[v], observed_value, &str);
374
375               /* The key */
376               tab_text (freq_table, 0, i + 1, TAB_LEFT, ds_cstr (&str));
377               ds_destroy (&str);
378
379
380               /* The observed N */
381               tab_float (freq_table, 1, i + 1, TAB_NONE,
382                          ff[i]->count, 8, 0);
383
384               if ( cst->n_expected > 0 )
385                 exp = cst->expected[i] * total_obs / total_expected ;
386               else
387                 exp = total_obs / (double) n_cells;
388
389               tab_float (freq_table, 2, i + 1, TAB_NONE,
390                          exp, 8, 2);
391
392               /* The residual */
393               tab_float (freq_table, 3, i + 1, TAB_NONE,
394                          ff[i]->count - exp, 8, 2);
395
396               xsq[v] += (ff[i]->count - exp) * (ff[i]->count - exp) / exp;
397             }
398
399           df[v] = n_cells - 1.0;
400
401           tab_float (freq_table, 1, i + 1, TAB_NONE,
402                      total_obs, 8, 0);
403
404           tab_submit (freq_table);
405
406           hsh_destroy (freq_hash);
407         }
408     }
409   else  /* ranged == true */
410     {
411       struct tab_table *freq_table = create_combo_frequency_table (cst);
412
413       n_cells = cst->hi - cst->lo + 1;
414
415       for ( v = 0 ; v < ost->n_vars ; ++v )
416         {
417           double total_obs = 0.0;
418           struct casereader *reader =
419             casereader_create_filter_missing (casereader_clone (input),
420                                               &ost->vars[v], 1, exclude,
421                                               NULL, NULL);
422           struct hsh_table *freq_hash =
423             create_freq_hash_with_range (dict, reader,
424                                          ost->vars[v], cst->lo, cst->hi);
425
426           struct freq **ff;
427
428           if (freq_hash == NULL)
429             continue;
430
431           ff = (struct freq **) hsh_sort (freq_hash);
432           assert ( n_cells == hsh_count (freq_hash));
433
434           for ( i = 0 ; i < hsh_count (freq_hash) ; ++i )
435             total_obs += ff[i]->count;
436
437           xsq[v] = 0.0;
438           for ( i = 0 ; i < hsh_count (freq_hash) ; ++i )
439             {
440               struct string str;
441               double exp;
442
443               const union value *observed_value = ff[i]->value;
444
445               ds_init_empty (&str);
446               var_append_value_name (ost->vars[v], observed_value, &str);
447               /* The key */
448               tab_text  (freq_table, v * 4 + 1, i + 2 , TAB_LEFT,
449                          ds_cstr (&str));
450               ds_destroy (&str);
451
452               /* The observed N */
453               tab_float (freq_table, v * 4 + 2, i + 2 , TAB_NONE,
454                          ff[i]->count, 8, 0);
455
456               if ( cst->n_expected > 0 )
457                 exp = cst->expected[i] * total_obs / total_expected ;
458               else
459                 exp = total_obs / (double) hsh_count (freq_hash);
460
461               /* The expected N */
462               tab_float (freq_table, v * 4 + 3, i + 2 , TAB_NONE,
463                          exp, 8, 2);
464
465               /* The residual */
466               tab_float (freq_table, v * 4 + 4, i + 2 , TAB_NONE,
467                          ff[i]->count - exp, 8, 2);
468
469               xsq[v] += (ff[i]->count - exp) * (ff[i]->count - exp) / exp;
470             }
471
472
473           tab_float (freq_table, v * 4 + 2, tab_nr (freq_table) - 1, TAB_NONE,
474                      total_obs, 8, 0);
475
476           df[v] = n_cells - 1.0;
477
478           hsh_destroy (freq_hash);
479         }
480
481       tab_submit (freq_table);
482     }
483   ok = !taint_has_tainted_successor (casereader_get_taint (input));
484   casereader_destroy (input);
485
486   if (ok)
487     {
488       struct tab_table *stats_table = create_stats_table (cst);
489
490       /* Populate the summary statistics table */
491       for ( v = 0 ; v < ost->n_vars ; ++v )
492         {
493           const struct variable *var = ost->vars[v];
494
495           tab_text (stats_table, 1 + v, 0, TAB_CENTER, var_get_name (var));
496
497           tab_float (stats_table, 1 + v, 1, TAB_NONE, xsq[v], 8,3);
498           tab_float (stats_table, 1 + v, 2, TAB_NONE, df[v], 8,0);
499
500           tab_float (stats_table, 1 + v, 3, TAB_NONE,
501                      gsl_cdf_chisq_Q (xsq[v], df[v]), 8,3);
502         }
503       tab_submit (stats_table);
504     }
505
506   free (xsq);
507   free (df);
508 }
509