1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2006, 2007, 2009, 2010, 2011 Free Software Foundation, Inc.
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.
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.
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/>. */
19 #include "language/stats/chisquare.h"
21 #include <gsl/gsl_cdf.h>
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"
43 #include "gl/xalloc.h"
46 #define _(msgid) gettext (msgid)
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. */
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)
60 struct freq **entries;
66 assert (var_is_numeric (var));
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 )
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));
80 for (; (c = casereader_read (input)) != NULL; case_unref (c))
82 double x = trunc (case_num (c, var));
83 if (x >= lo && x <= hi)
86 struct freq *fr = entries[ofs];
87 fr->count += dict_get_case_weight (dict, c, &warn);
93 return casereader_destroy (input);
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
102 create_freq_hash (const struct dictionary *dict,
103 struct casereader *input,
104 const struct variable *var,
105 struct hmap *freq_hash)
107 int width = var_get_width (var);
111 for (; (c = casereader_read (input)) != NULL; case_unref (c))
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);
118 f = freq_hmap_search (freq_hash, value, width, hash);
120 f = freq_hmap_insert (freq_hash, value, width, hash);
125 return casereader_destroy (input);
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)
136 const struct one_sample_test *ost = (const struct one_sample_test*)test;
138 struct tab_table *table ;
139 const struct variable *var = ost->vars[v];
141 hmap_init (freq_hash);
142 if (!create_freq_hash (dict, input, var, freq_hash))
144 freq_hmap_destroy (freq_hash, var_get_width (var));
148 n_cells = hmap_count (freq_hash);
150 if ( test->n_expected > 0 && n_cells != test->n_expected )
152 msg(ME, _("CHISQUARE test specified %d expected values, but"
153 " %d distinct values were encountered in variable %s."),
154 test->n_expected, n_cells,
160 table = tab_create(4, n_cells + 2);
162 tab_title (table, "%s", var_to_string(var));
163 tab_text (table, 1, 0, TAB_LEFT, _("Observed N"));
164 tab_text (table, 2, 0, TAB_LEFT, _("Expected N"));
165 tab_text (table, 3, 0, TAB_LEFT, _("Residual"));
167 tab_headers (table, 1, 0, 1, 0);
169 tab_box (table, TAL_1, TAL_1, -1, -1,
170 0, 0, tab_nc (table) - 1, tab_nr(table) - 1 );
172 tab_hline (table, TAL_1, 0, tab_nc(table) - 1, 1);
174 tab_vline (table, TAL_2, 1, 0, tab_nr(table) - 1);
175 for ( i = 2 ; i < 4 ; ++i )
176 tab_vline (table, TAL_1, i, 0, tab_nr(table) - 1);
179 tab_text (table, 0, tab_nr (table) - 1, TAB_LEFT, _("Total"));
185 static struct tab_table *
186 create_combo_frequency_table (const struct chisquare_test *test)
189 const struct one_sample_test *ost = (const struct one_sample_test*)test;
191 struct tab_table *table ;
193 int n_cells = test->hi - test->lo + 1;
195 table = tab_create(1 + ost->n_vars * 4, n_cells + 3);
197 tab_title (table, _("Frequencies"));
198 for ( i = 0 ; i < ost->n_vars ; ++i )
200 const struct variable *var = ost->vars[i];
201 tab_text (table, i * 4 + 1, 1, TAB_LEFT, _("Category"));
202 tab_text (table, i * 4 + 2, 1, TAB_LEFT, _("Observed N"));
203 tab_text (table, i * 4 + 3, 1, TAB_LEFT, _("Expected N"));
204 tab_text (table, i * 4 + 4, 1, TAB_LEFT, _("Residual"));
206 tab_vline (table, TAL_2, i * 4 + 1,
207 0, tab_nr (table) - 1);
209 tab_vline (table, TAL_1, i * 4 + 2,
210 0, tab_nr (table) - 1);
212 tab_vline (table, TAL_1, i * 4 + 3,
213 1, tab_nr (table) - 1);
215 tab_vline (table, TAL_1, i * 4 + 4,
216 1, tab_nr (table) - 1);
219 tab_joint_text (table,
223 var_to_string (var));
226 for ( i = test->lo ; i <= test->hi ; ++i )
227 tab_fixed (table, 0, 2 + i - test->lo,
228 TAB_LEFT, 1 + i - test->lo, 8, 0);
230 tab_headers (table, 1, 0, 2, 0);
232 tab_box (table, TAL_1, TAL_1, -1, -1,
233 0, 0, tab_nc (table) - 1, tab_nr(table) - 1 );
235 tab_hline (table, TAL_1, 1, tab_nc(table) - 1, 1);
236 tab_hline (table, TAL_1, 0, tab_nc(table) - 1, 2);
238 tab_text (table, 0, tab_nr (table) - 1, TAB_LEFT, _("Total"));
244 static struct tab_table *
245 create_stats_table (const struct chisquare_test *test)
247 const struct one_sample_test *ost = (const struct one_sample_test*) test;
249 struct tab_table *table;
250 table = tab_create (1 + ost->n_vars, 4);
251 tab_title (table, _("Test Statistics"));
252 tab_headers (table, 1, 0, 1, 0);
254 tab_box (table, TAL_1, TAL_1, -1, -1,
255 0, 0, tab_nc(table) - 1, tab_nr(table) - 1 );
257 tab_box (table, -1, -1, -1, TAL_1,
258 1, 0, tab_nc(table) - 1, tab_nr(table) - 1 );
261 tab_vline (table, TAL_2, 1, 0, tab_nr (table) - 1);
262 tab_hline (table, TAL_1, 0, tab_nc (table) - 1, 1);
265 tab_text (table, 0, 1, TAB_LEFT, _("Chi-Square"));
266 tab_text (table, 0, 2, TAB_LEFT, _("df"));
267 tab_text (table, 0, 3, TAB_LEFT, _("Asymp. Sig."));
274 chisquare_execute (const struct dataset *ds,
275 struct casereader *input,
276 enum mv_class exclude,
277 const struct npar_test *test,
281 const struct dictionary *dict = dataset_dict (ds);
283 struct chisquare_test *cst = UP_CAST (test, struct chisquare_test,
285 struct one_sample_test *ost = &cst->parent;
287 double total_expected = 0.0;
288 const struct variable *wvar = dict_get_weight (dict);
289 const struct fmt_spec *wfmt = wvar ?
290 var_get_print_format (wvar) : & F_8_0;
292 double *df = xzalloc (sizeof (*df) * ost->n_vars);
293 double *xsq = xzalloc (sizeof (*df) * ost->n_vars);
296 for ( i = 0 ; i < cst->n_expected ; ++i )
297 total_expected += cst->expected[i];
299 if ( cst->ranged == false )
301 for ( v = 0 ; v < ost->n_vars ; ++v )
303 const struct variable *var = ost->vars[v];
304 double total_obs = 0.0;
305 struct hmap freq_hash;
306 struct casereader *reader =
307 casereader_create_filter_missing (casereader_clone (input),
310 struct tab_table *freq_table =
311 create_variable_frequency_table(dict, reader, cst, v, &freq_hash);
315 if ( NULL == freq_table )
317 ff = freq_hmap_sort (&freq_hash, var_get_width (var));
319 n_cells = hmap_count (&freq_hash);
321 for ( i = 0 ; i < n_cells ; ++i )
322 total_obs += ff[i]->count;
325 for ( i = 0 ; i < n_cells ; ++i )
329 const union value *observed_value = &ff[i]->value;
331 ds_init_empty (&str);
332 var_append_value_name (var, observed_value, &str);
335 tab_text (freq_table, 0, i + 1, TAB_LEFT, ds_cstr (&str));
340 tab_double (freq_table, 1, i + 1, TAB_NONE,
343 if ( cst->n_expected > 0 )
344 exp = cst->expected[i] * total_obs / total_expected ;
346 exp = total_obs / (double) n_cells;
348 tab_double (freq_table, 2, i + 1, TAB_NONE,
352 tab_double (freq_table, 3, i + 1, TAB_NONE,
353 ff[i]->count - exp, NULL);
355 xsq[v] += (ff[i]->count - exp) * (ff[i]->count - exp) / exp;
358 df[v] = n_cells - 1.0;
360 tab_double (freq_table, 1, i + 1, TAB_NONE,
363 tab_submit (freq_table);
365 freq_hmap_destroy (&freq_hash, var_get_width (var));
369 else /* ranged == true */
371 struct tab_table *freq_table = create_combo_frequency_table (cst);
373 n_cells = cst->hi - cst->lo + 1;
375 for ( v = 0 ; v < ost->n_vars ; ++v )
377 const struct variable *var = ost->vars[v];
378 double total_obs = 0.0;
379 struct casereader *reader =
380 casereader_create_filter_missing (casereader_clone (input),
383 struct hmap freq_hash;
386 hmap_init (&freq_hash);
387 if (!create_freq_hash_with_range (dict, reader, var,
388 cst->lo, cst->hi, &freq_hash))
390 freq_hmap_destroy (&freq_hash, var_get_width (var));
394 ff = freq_hmap_sort (&freq_hash, var_get_width (var));
396 for ( i = 0 ; i < hmap_count (&freq_hash) ; ++i )
397 total_obs += ff[i]->count;
400 for ( i = 0 ; i < hmap_count (&freq_hash) ; ++i )
405 const union value *observed_value = &ff[i]->value;
407 ds_init_empty (&str);
408 var_append_value_name (ost->vars[v], observed_value, &str);
410 tab_text (freq_table, v * 4 + 1, i + 2 , TAB_LEFT,
415 tab_double (freq_table, v * 4 + 2, i + 2 , TAB_NONE,
418 if ( cst->n_expected > 0 )
419 exp = cst->expected[i] * total_obs / total_expected ;
421 exp = total_obs / (double) hmap_count (&freq_hash);
424 tab_double (freq_table, v * 4 + 3, i + 2 , TAB_NONE,
428 tab_double (freq_table, v * 4 + 4, i + 2 , TAB_NONE,
429 ff[i]->count - exp, NULL);
431 xsq[v] += (ff[i]->count - exp) * (ff[i]->count - exp) / exp;
435 tab_double (freq_table, v * 4 + 2, tab_nr (freq_table) - 1, TAB_NONE,
438 df[v] = n_cells - 1.0;
440 freq_hmap_destroy (&freq_hash, var_get_width (var));
444 tab_submit (freq_table);
446 ok = !taint_has_tainted_successor (casereader_get_taint (input));
447 casereader_destroy (input);
451 struct tab_table *stats_table = create_stats_table (cst);
453 /* Populate the summary statistics table */
454 for ( v = 0 ; v < ost->n_vars ; ++v )
456 const struct variable *var = ost->vars[v];
458 tab_text (stats_table, 1 + v, 0, TAB_CENTER, var_get_name (var));
460 tab_double (stats_table, 1 + v, 1, TAB_NONE, xsq[v], NULL);
461 tab_fixed (stats_table, 1 + v, 2, TAB_NONE, df[v], 8, 0);
463 tab_double (stats_table, 1 + v, 3, TAB_NONE,
464 gsl_cdf_chisq_Q (xsq[v], df[v]), NULL);
466 tab_submit (stats_table);