1 /* PSPP - computes sample statistics.
2 Copyright (C) 2006 Free Software Foundation, Inc.
3 Written by John Darrington <john@darrington.wattle.id.au>
5 This program is free software; you can redistribute it and/or
6 modify it under the terms of the GNU General Public License as
7 published by the Free Software Foundation; either version 2 of the
8 License, or (at your option) any later version.
10 This program is distributed in the hope that it will be useful, but
11 WITHOUT ANY WARRANTY; without even the implied warranty of
12 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 General Public License for more details.
15 You should have received a copy of the GNU General Public License
16 along with this program; if not, write to the Free Software
17 Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
21 #include <libpspp/compiler.h>
22 #include <output/table.h>
23 #include <libpspp/alloc.h>
25 #include <data/case.h>
26 #include <data/casefile.h>
27 #include <data/dictionary.h>
28 #include <data/procedure.h>
29 #include <data/variable.h>
30 #include <data/value.h>
31 #include <data/value-labels.h>
32 #include <data/casefilter.h>
34 #include <libpspp/message.h>
35 #include <libpspp/assertion.h>
41 #define _(msgid) gettext (msgid)
43 #include <libpspp/misc.h>
45 #include <gsl/gsl_cdf.h>
46 #include <gsl/gsl_randist.h>
47 #include <gsl-extras/gsl-extras.h>
51 #include <libpspp/hash.h>
53 static double calculate_binomial_internal (double n1, double n2,
58 swap (double *i1, double *i2)
66 calculate_binomial (double n1, double n2, double p)
68 const double n = n1 + n2;
69 const bool test_reversed = (n1 / n > p ) ;
76 return calculate_binomial_internal (n1, n2, p);
80 calculate_binomial_internal (double n1, double n2, double p)
82 /* SPSS Statistical Algorithms has completely different and WRONG
85 double sig1tailed = gslextras_cdf_binomial_P (n1, n1 + n2, p);
88 return sig1tailed > 0.5 ? 1.0 :sig1tailed * 2.0;
94 do_binomial (const struct dictionary *dict,
95 const struct casefile *cf,
96 const struct binomial_test *bst,
99 const struct casefilter *filter
104 const struct one_sample_test *ost = (const struct one_sample_test *) bst;
106 struct casereader *r = casefile_get_reader (cf, NULL);
108 while (casereader_read(r, &c))
112 dict_get_case_weight (dict, &c, &warn);
114 for (v = 0 ; v < ost->n_vars ; ++v )
116 const struct variable *var = ost->vars[v];
117 const union value *value = case_data (&c, var);
119 if ( casefilter_variable_missing (filter, &c, var))
122 if ( NULL == cat1[v].value )
124 cat1[v].value = value_dup (value, var_get_width (var));
127 else if ( 0 == compare_values (cat1[v].value, value,
128 var_get_width (var)))
130 else if ( NULL == cat2[v].value )
132 cat2[v].value = value_dup (value, var_get_width (var));
135 else if ( 0 == compare_values (cat2[v].value, value,
136 var_get_width (var)))
138 else if ( bst->category1 == SYSMIS)
139 msg (ME, _("Variable %s is not dichotomous"), var_get_name (var));
144 casereader_destroy (r);
150 binomial_execute (const struct dataset *ds,
151 const struct casefile *cf,
152 struct casefilter *filter,
153 const struct npar_test *test)
156 const struct binomial_test *bst = (const struct binomial_test *) test;
157 const struct one_sample_test *ost = (const struct one_sample_test*) test;
159 struct freq *cat1 = xzalloc (sizeof (*cat1) * ost->n_vars);
160 struct freq *cat2 = xzalloc (sizeof (*cat1) * ost->n_vars);
161 struct tab_table *table ;
163 assert ((bst->category1 == SYSMIS) == (bst->category2 == SYSMIS) );
165 if ( bst->category1 != SYSMIS )
168 v.f = bst->category1;
169 cat1->value = value_dup (&v, 0);
172 if ( bst->category2 != SYSMIS )
175 v.f = bst->category2;
176 cat2->value = value_dup (&v, 0);
179 do_binomial (dataset_dict(ds), cf, bst, cat1, cat2, filter);
181 table = tab_create (7, ost->n_vars * 3 + 1, 0);
183 tab_dim (table, tab_natural_dimensions);
185 tab_title (table, _("Binomial Test"));
187 tab_headers (table, 2, 0, 1, 0);
189 tab_box (table, TAL_1, TAL_1, -1, TAL_1,
190 0, 0, table->nc - 1, tab_nr(table) - 1 );
192 for (v = 0 ; v < ost->n_vars; ++v)
195 const struct variable *var = ost->vars[v];
196 tab_hline (table, TAL_1, 0, tab_nc (table) -1, 1 + v * 3);
199 tab_text (table, 0, 1 + v * 3, TAB_LEFT,
200 var_to_string (var));
202 tab_text (table, 1, 1 + v * 3, TAB_LEFT,
205 tab_text (table, 1, 2 + v * 3, TAB_LEFT,
208 tab_text (table, 1, 3 + v * 3, TAB_LEFT,
212 tab_float (table, 5, 1 + v * 3, TAB_NONE, bst->p, 8, 3);
214 /* Category labels */
215 tab_text (table, 2, 1 + v * 3, TAB_NONE,
216 var_get_value_name (var, cat1[v].value));
218 tab_text (table, 2, 2 + v * 3, TAB_NONE,
219 var_get_value_name (var, cat2[v].value));
222 tab_float (table, 3, 1 + v * 3, TAB_NONE,
223 cat1[v].count, 8, 0);
225 tab_float (table, 3, 2 + v * 3, TAB_NONE,
226 cat2[v].count, 8, 0);
228 n_total = cat1[v].count + cat2[v].count;
231 tab_float (table, 3, 3 + v * 3, TAB_NONE,
234 /* Observed Proportions */
236 tab_float (table, 4, 1 + v * 3, TAB_NONE,
237 cat1[v].count / n_total, 8, 3);
239 tab_float (table, 4, 2 + v * 3, TAB_NONE,
240 cat2[v].count / n_total, 8, 3);
242 tab_float (table, 4, 3 + v * 3, TAB_NONE,
243 (cat1[v].count + cat2[v].count) / n_total, 8, 2);
247 sig = calculate_binomial (cat1[v].count, cat2[v].count,
250 tab_float (table, 6, 1 + v * 3, TAB_NONE,
254 tab_text (table, 2, 0, TAB_CENTER, _("Category"));
255 tab_text (table, 3, 0, TAB_CENTER, _("N"));
256 tab_text (table, 4, 0, TAB_CENTER, _("Observed Prop."));
257 tab_text (table, 5, 0, TAB_CENTER, _("Test Prop."));
259 tab_text (table, 6, 0, TAB_CENTER | TAT_PRINTF,
260 _("Exact Sig. (%d-tailed)"),
261 bst->p == 0.5 ? 2: 1);
263 tab_vline (table, TAL_2, 2, 0, tab_nr (table) -1);