Update all #include directives to the currently preferred style.
[pspp-builds.git] / src / language / stats / cochran.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 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/cochran.h"
20
21 #include <gsl/gsl_cdf.h>
22 #include <stdbool.h>
23
24 #include "data/casereader.h"
25 #include "data/dictionary.h"
26 #include "data/format.h"
27 #include "data/procedure.h"
28 #include "data/val-type.h"
29 #include "data/variable.h"
30 #include "language/stats/npar.h"
31 #include "libpspp/cast.h"
32 #include "libpspp/message.h"
33 #include "libpspp/misc.h"
34 #include "output/tab.h"
35
36 #include "gettext.h"
37 #define _(msgid) gettext (msgid)
38
39 struct cochran
40 {
41   double success;
42   double failure;
43
44   double *hits;
45   double *misses;
46
47   const struct dictionary *dict;
48   double cc;
49   double df;
50   double q;
51 };
52
53 static void show_freqs_box (const struct one_sample_test *ost, const struct cochran *ch);
54 static void show_sig_box (const struct cochran *ch);
55
56 void
57 cochran_execute (const struct dataset *ds,
58               struct casereader *input,
59               enum mv_class exclude,
60               const struct npar_test *test,
61               bool exact UNUSED, double timer UNUSED)
62 {
63   struct one_sample_test *ct = UP_CAST (test, struct one_sample_test, parent);
64   int v;
65   struct cochran ch;
66   const struct dictionary *dict = dataset_dict (ds);
67   const struct variable *weight = dict_get_weight (dict);
68
69   struct ccase *c;
70   double rowsq = 0;
71   ch.cc = 0.0;
72   ch.dict = dict;
73   ch.success = SYSMIS;
74   ch.failure = SYSMIS;
75   ch.hits = xcalloc (ct->n_vars, sizeof *ch.hits);
76   ch.misses = xcalloc (ct->n_vars, sizeof *ch.misses);
77
78   for (; (c = casereader_read (input)); case_unref (c))
79     {
80       double case_hits = 0.0;
81       const double w = weight ? case_data (c, weight)->f: 1.0;
82       for (v = 0; v < ct->n_vars; ++v)
83         {
84           const struct variable *var = ct->vars[v];
85           const union value *val = case_data (c, var);
86
87           if ( var_is_value_missing (var, val, exclude))
88             continue;
89
90           if ( ch.success == SYSMIS)
91             {
92               ch.success = val->f;
93             }
94           else if (ch.failure == SYSMIS && val->f != ch.success)
95             {
96               ch.failure = val->f;
97             }
98           if ( ch.success == val->f)
99             {
100               ch.hits[v] += w;
101               case_hits += w;
102             }
103           else if ( ch.failure == val->f)
104             {
105               ch.misses[v] += w;
106             }
107           else
108             {
109               msg (MW, _("More than two values encountered.  Cochran Q test will not be run."));
110               goto finish;
111             }
112         }
113       ch.cc += w;
114       rowsq += pow2 (case_hits);
115     }
116   casereader_destroy (input);
117   
118   {
119     double c_l = 0;
120     double c_l2 = 0;
121     for (v = 0; v < ct->n_vars; ++v)
122       {
123         c_l += ch.hits[v];
124         c_l2 += pow2 (ch.hits[v]);
125       }
126
127     ch.q = ct->n_vars * c_l2;
128     ch.q -= pow2 (c_l);
129     ch.q *= ct->n_vars - 1;
130
131     ch.q /= ct->n_vars * c_l - rowsq;
132   
133     ch.df = ct->n_vars - 1;
134   }
135
136   show_freqs_box (ct, &ch);
137   show_sig_box (&ch);
138
139  finish:
140
141   free (ch.hits);
142   free (ch.misses);
143 }
144
145 static void
146 show_freqs_box (const struct one_sample_test *ost, const struct cochran *ct)
147 {
148   int i;
149   const struct variable *weight = dict_get_weight (ct->dict);
150   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
151
152   const int row_headers = 1;
153   const int column_headers = 2;
154   struct tab_table *table =
155     tab_create (row_headers + 2, column_headers + ost->n_vars);
156
157   tab_headers (table, row_headers, 0, column_headers, 0);
158
159   tab_title (table, _("Frequencies"));
160
161   /* Vertical lines inside the box */
162   tab_box (table, 1, 0, -1, TAL_1,
163            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
164
165   /* Box around the table */
166   tab_box (table, TAL_2, TAL_2, -1, -1,
167            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
168
169   tab_joint_text (table, 1, 0, 2, 0,
170                   TAT_TITLE | TAB_CENTER, _("Value"));
171
172   tab_text_format (table, 1, 1, 0, _("Success (%g)"), ct->success);
173   tab_text_format (table, 2, 1, 0, _("Failure (%g)"), ct->failure);
174
175   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, column_headers);
176   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
177
178   for (i = 0 ; i < ost->n_vars ; ++i)
179     {
180       tab_text (table, 0, column_headers + i,
181                 TAB_LEFT, var_to_string (ost->vars[i]));
182
183       tab_double (table, 1, column_headers + i, 0,
184                   ct->hits[i], wfmt);
185
186       tab_double (table, 2, column_headers + i, 0,
187                   ct->misses[i], wfmt);
188     }
189
190   tab_submit (table);
191 }
192
193
194
195 static void
196 show_sig_box (const struct cochran *ch)
197 {
198   const struct variable *weight = dict_get_weight (ch->dict);
199   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
200
201   const int row_headers = 1;
202   const int column_headers = 0;
203   struct tab_table *table =
204     tab_create (row_headers + 1, column_headers + 4);
205
206   tab_headers (table, row_headers, 0, column_headers, 0);
207
208   tab_title (table, _("Test Statistics"));
209
210   tab_text (table,  0, column_headers,
211             TAT_TITLE | TAB_LEFT , _("N"));
212
213   tab_text (table,  0, 1 + column_headers,
214             TAT_TITLE | TAB_LEFT , _("Cochran's Q"));
215
216   tab_text (table,  0, 2 + column_headers,
217             TAT_TITLE | TAB_LEFT, _("df"));
218
219   tab_text (table,  0, 3 + column_headers,
220             TAT_TITLE | TAB_LEFT, _("Asymp. Sig."));
221
222   /* Box around the table */
223   tab_box (table, TAL_2, TAL_2, -1, -1,
224            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
225
226   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
227   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
228
229   tab_double (table, 1, column_headers, 
230               0, ch->cc, wfmt);
231
232   tab_double (table, 1, column_headers + 1, 
233               0, ch->q, 0);
234
235   tab_double (table, 1, column_headers + 2, 
236               0, ch->df, &F_8_0);
237
238   tab_double (table, 1, column_headers + 3, 
239               0, gsl_cdf_chisq_Q (ch->q, ch->df), 
240               0);
241
242   tab_submit (table);
243 }