GLM: Rewrite interactions module and update glm.c to use it in some places
[pspp] / src / language / stats / sign.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2009, 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/sign.h"
20
21 #include <gsl/gsl_cdf.h>
22 #include <gsl/gsl_randist.h>
23
24 #include "data/casereader.h"
25 #include "data/dataset.h"
26 #include "data/dictionary.h"
27 #include "data/format.h"
28 #include "data/missing-values.h"
29 #include "data/variable.h"
30 #include "language/stats/npar.h"
31 #include "libpspp/str.h"
32 #include "output/tab.h"
33
34 #include "gl/minmax.h"
35 #include "gl/xalloc.h"
36
37 #include "gettext.h"
38 #define _(msgid) gettext (msgid)
39
40 struct sign_test_params
41 {
42   double pos;
43   double ties;
44   double neg;
45
46   double one_tailed_sig;
47   double point_prob;
48 };
49
50
51 static void
52 output_frequency_table (const struct two_sample_test *t2s,
53                         const struct sign_test_params *param,
54                         const struct dictionary *dict)
55 {
56   int i;
57   struct tab_table *table = tab_create (3, 1 + 4 * t2s->n_pairs);
58
59   const struct variable *wv = dict_get_weight (dict);
60   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : & F_8_0;
61
62   tab_title (table, _("Frequencies"));
63
64   tab_headers (table, 2, 0, 1, 0);
65
66   /* Vertical lines inside the box */
67   tab_box (table, 0, 0, -1, TAL_1,
68            1, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
69
70   /* Box around entire table */
71   tab_box (table, TAL_2, TAL_2, -1, -1,
72            0, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
73
74   tab_text (table,  2, 0,  TAB_CENTER, _("N"));
75
76   for (i = 0 ; i < t2s->n_pairs; ++i)
77     {
78       variable_pair *vp = &t2s->pairs[i];
79
80       struct string pair_name;
81       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
82       ds_put_cstr (&pair_name, " - ");
83       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
84
85
86       tab_text (table, 0, 1 + i * 4, TAB_LEFT, ds_cstr (&pair_name));
87
88       ds_destroy (&pair_name);
89
90       tab_hline (table, TAL_1, 0, tab_nc (table) - 1, 1 + i * 4);
91
92       tab_text (table,  1, 1 + i * 4,  TAB_LEFT, _("Negative Differences"));
93       tab_text (table,  1, 2 + i * 4,  TAB_LEFT, _("Positive Differences"));
94       tab_text (table,  1, 3 + i * 4,  TAB_LEFT, _("Ties"));
95       tab_text (table,  1, 4 + i * 4,  TAB_LEFT, _("Total"));
96
97       tab_double (table, 2, 1 + i * 4, TAB_RIGHT, param[i].neg, wfmt);
98       tab_double (table, 2, 2 + i * 4, TAB_RIGHT, param[i].pos, wfmt);
99       tab_double (table, 2, 3 + i * 4, TAB_RIGHT, param[i].ties, wfmt);
100       tab_double (table, 2, 4 + i * 4, TAB_RIGHT,
101                  param[i].ties + param[i].neg + param[i].pos, wfmt);
102     }
103
104   tab_submit (table);
105 }
106
107 static void
108 output_statistics_table (const struct two_sample_test *t2s,
109                          const struct sign_test_params *param)
110 {
111   int i;
112   struct tab_table *table = tab_create (1 + t2s->n_pairs, 4);
113
114   tab_title (table, _("Test Statistics"));
115
116   tab_headers (table, 0, 1,  0, 1);
117
118   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, 1);
119   tab_vline (table, TAL_2, 1, 0, tab_nr (table) - 1);
120
121
122   /* Vertical lines inside the box */
123   tab_box (table, -1, -1, -1, TAL_1,
124            0, 0,
125            tab_nc (table) - 1, tab_nr (table) - 1);
126
127   /* Box around entire table */
128   tab_box (table, TAL_2, TAL_2, -1, -1,
129            0, 0, tab_nc (table) - 1,
130            tab_nr (table) - 1);
131
132   tab_text (table,  0, 1, TAT_TITLE | TAB_LEFT,
133             _("Exact Sig. (2-tailed)"));
134
135   tab_text (table,  0, 2, TAT_TITLE | TAB_LEFT,
136             _("Exact Sig. (1-tailed)"));
137
138   tab_text (table,  0, 3, TAT_TITLE | TAB_LEFT,
139             _("Point Probability"));
140
141   for (i = 0 ; i < t2s->n_pairs; ++i)
142     {
143       variable_pair *vp = &t2s->pairs[i];
144
145       struct string pair_name;
146       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
147       ds_put_cstr (&pair_name, " - ");
148       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
149
150       tab_text (table,  1 + i, 0, TAB_LEFT, ds_cstr (&pair_name));
151       ds_destroy (&pair_name);
152
153       tab_double (table, 1 + i, 1, TAB_RIGHT,
154                   param[i].one_tailed_sig * 2, NULL);
155
156       tab_double (table, 1 + i, 2, TAB_RIGHT, param[i].one_tailed_sig, NULL);
157       tab_double (table, 1 + i, 3, TAB_RIGHT, param[i].point_prob, NULL);
158     }
159
160   tab_submit (table);
161 }
162
163 void
164 sign_execute (const struct dataset *ds,
165                   struct casereader *input,
166                   enum mv_class exclude,
167                   const struct npar_test *test,
168                   bool exact UNUSED,
169                   double timer UNUSED)
170 {
171   int i;
172   bool warn = true;
173   const struct dictionary *dict = dataset_dict (ds);
174   const struct two_sample_test *t2s = UP_CAST (test, const struct two_sample_test, parent);
175   struct ccase *c;
176
177   struct sign_test_params *stp = xcalloc (sizeof *stp, t2s->n_pairs);
178
179   struct casereader *r = input;
180
181   for (; (c = casereader_read (r)) != NULL; case_unref (c))
182     {
183       const double weight = dict_get_case_weight (dict, c, &warn);
184
185       for (i = 0 ; i < t2s->n_pairs; ++i )
186         {
187           variable_pair *vp = &t2s->pairs[i];
188           const union value *value0 = case_data (c, (*vp)[0]);
189           const union value *value1 = case_data (c, (*vp)[1]);
190           const double diff = value0->f - value1->f;
191
192           if (var_is_value_missing ((*vp)[0], value0, exclude))
193             continue;
194
195           if (var_is_value_missing ((*vp)[1], value1, exclude))
196             continue;
197
198           if ( diff > 0)
199             stp[i].pos += weight;
200           else if (diff < 0)
201             stp[i].neg += weight;
202           else
203             stp[i].ties += weight;
204         }
205     }
206
207   casereader_destroy (r);
208
209   for (i = 0 ; i < t2s->n_pairs; ++i )
210     {
211       int r = MIN (stp[i].pos, stp[i].neg);
212       stp[i].one_tailed_sig = gsl_cdf_binomial_P (r,
213                                                   0.5,
214                                                   stp[i].pos + stp[i].neg);
215
216       stp[i].point_prob = gsl_ran_binomial_pdf (r, 0.5,
217                                                 stp[i].pos + stp[i].neg);
218     }
219
220   output_frequency_table (t2s, stp, dict);
221
222   output_statistics_table (t2s, stp);
223
224   free (stp);
225 }