expressions: Implement the STRUNC function.
[pspp] / src / language / stats / mann-whitney.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/mann-whitney.h"
20
21 #include <gsl/gsl_cdf.h>
22
23 #include "data/case.h"
24 #include "data/casereader.h"
25 #include "data/dataset.h"
26 #include "data/dictionary.h"
27 #include "data/format.h"
28 #include "data/variable.h"
29 #include "libpspp/cast.h"
30 #include "libpspp/misc.h"
31 #include "math/sort.h"
32 #include "output/tab.h"
33
34 /* Calculates the adjustment necessary for tie compensation */
35 static void
36 distinct_callback (double v UNUSED, casenumber t, double w UNUSED, void *aux)
37 {
38   double *tiebreaker = aux;
39
40   *tiebreaker += (pow3 (t) - t) / 12.0;
41 }
42
43 struct mw
44 {
45   double rank_sum[2];
46   double n[2];
47
48   double u;  /* The Mann-Whitney U statistic */
49   double w;  /* The Wilcoxon Rank Sum W statistic */
50   double z;  
51 };
52
53 static void show_ranks_box (const struct n_sample_test *nst, const struct mw *mw);
54 static void show_statistics_box (const struct n_sample_test *nst, const struct mw *mw, bool exact);
55
56
57
58 static bool
59 belongs_to_test (const struct ccase *c, void *aux)
60 {
61   const struct n_sample_test *nst = aux;
62
63   const union value *group = case_data (c, nst->indep_var);
64   const size_t group_var_width = var_get_width (nst->indep_var);
65
66   if ( value_equal (group, &nst->val1, group_var_width))
67     return true;
68
69   if ( value_equal (group, &nst->val2, group_var_width))
70     return true;
71
72   return false;
73 }
74
75                                          
76
77 void
78 mann_whitney_execute (const struct dataset *ds,
79                       struct casereader *input,
80                       enum mv_class exclude,
81                       const struct npar_test *test,
82                       bool exact,
83                       double timer UNUSED)
84 {
85   int i;
86   const struct dictionary *dict = dataset_dict (ds);
87   const struct n_sample_test *nst = UP_CAST (test, const struct n_sample_test, parent);
88
89   const struct caseproto *proto = casereader_get_proto (input);
90   size_t rank_idx = caseproto_get_n_widths (proto);
91
92   struct mw *mw = xcalloc (nst->n_vars, sizeof *mw);
93
94   for (i = 0; i < nst->n_vars; ++i)
95     {
96       double tiebreaker = 0.0;
97       bool warn = true;
98       enum rank_error rerr = 0;
99       struct casereader *rr;
100       struct ccase *c;
101       const struct variable *var = nst->vars[i];
102
103       struct casereader *reader = 
104         casereader_create_filter_func (casereader_clone (input),
105                                        belongs_to_test,
106                                        NULL,
107                                        CONST_CAST (struct n_sample_test *, nst),
108                                        NULL);
109
110       
111       reader = sort_execute_1var (reader, var);
112
113       rr = casereader_create_append_rank (reader, var,
114                                           dict_get_weight (dict),
115                                           &rerr,
116                                           distinct_callback, &tiebreaker);
117
118       for (; (c = casereader_read (rr)); case_unref (c))
119         {
120           const union value *val = case_data (c, var);
121           const union value *group = case_data (c, nst->indep_var);
122           const size_t group_var_width = var_get_width (nst->indep_var);
123           const double rank = case_data_idx (c, rank_idx)->f;
124
125           if ( var_is_value_missing (var, val, exclude))
126             continue;
127
128           if ( value_equal (group, &nst->val1, group_var_width))
129             {
130               mw[i].rank_sum[0] += rank;
131               mw[i].n[0] += dict_get_case_weight (dict, c, &warn);
132             }
133           else if ( value_equal (group, &nst->val2, group_var_width))
134             {
135               mw[i].rank_sum[1] += rank;
136               mw[i].n[1] += dict_get_case_weight (dict, c, &warn);
137             }
138         }
139       casereader_destroy (rr);
140
141       {
142         double n;
143         double denominator;
144         struct mw *mwv = &mw[i];
145
146         mwv->u = mwv->n[0] * mwv->n[1] ;
147         mwv->u += mwv->n[0] * (mwv->n[0] + 1) / 2.0;
148         mwv->u -= mwv->rank_sum[0];
149
150         mwv->w = mwv->rank_sum[1];
151         if ( mwv->u > mwv->n[0] * mwv->n[1] / 2.0)
152           {
153             mwv->u =  mwv->n[0] * mwv->n[1] - mwv->u;
154             mwv->w = mwv->rank_sum[0];
155           }
156         mwv->z = mwv->u - mwv->n[0] * mwv->n[1] / 2.0;
157         n = mwv->n[0] + mwv->n[1];
158         denominator = pow3(n) - n;
159         denominator /= 12;
160         denominator -= tiebreaker;
161         denominator *= mwv->n[0] * mwv->n[1];
162         denominator /= n * (n - 1);
163       
164         mwv->z /= sqrt (denominator);
165       }
166     }
167   casereader_destroy (input);
168
169   show_ranks_box (nst, mw);
170   show_statistics_box (nst, mw, exact);
171
172   free (mw);
173 }
174
175 \f
176
177 #include "gettext.h"
178 #define _(msgid) gettext (msgid)
179
180 static void
181 show_ranks_box (const struct n_sample_test *nst, const struct mw *mwv)
182 {
183   int i;
184   const int row_headers = 1;
185   const int column_headers = 2;
186   struct tab_table *table =
187     tab_create (row_headers + 7, column_headers + nst->n_vars);
188
189   struct string g1str, g2str;;
190   ds_init_empty (&g1str);
191   var_append_value_name (nst->indep_var, &nst->val1, &g1str);
192
193   ds_init_empty (&g2str);
194   var_append_value_name (nst->indep_var, &nst->val2, &g2str);
195
196   tab_headers (table, row_headers, 0, column_headers, 0);
197
198   tab_title (table, _("Ranks"));
199
200   /* Vertical lines inside the box */
201   tab_box (table, 1, 0, -1, TAL_1,
202            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
203
204   /* Box around the table */
205   tab_box (table, TAL_2, TAL_2, -1, -1,
206            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
207
208   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
209   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
210
211   tab_hline (table, TAL_1, row_headers, tab_nc (table) -1, 1);
212
213   tab_text (table, 1, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g1str));
214   tab_text (table, 2, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g2str));
215   tab_text (table, 3, 1, TAT_TITLE | TAB_CENTER, _("Total"));
216   tab_joint_text (table, 1, 0, 3, 0,
217                   TAT_TITLE | TAB_CENTER, _("N"));
218   tab_vline (table, TAL_2, 4, 0, tab_nr (table) - 1);
219
220   tab_text (table, 4, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g1str));
221   tab_text (table, 5, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g2str));
222   tab_joint_text (table, 4, 0, 5, 0,
223                   TAT_TITLE | TAB_CENTER, _("Mean Rank"));
224   tab_vline (table, TAL_2, 6, 0, tab_nr (table) - 1);
225
226   tab_text (table, 6, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g1str));
227   tab_text (table, 7, 1, TAT_TITLE | TAB_CENTER, ds_cstr (&g2str));
228   tab_joint_text (table, 6, 0, 7, 0,
229                   TAT_TITLE | TAB_CENTER, _("Sum of Ranks"));
230
231   ds_destroy (&g1str);
232   ds_destroy (&g2str);
233
234   for (i = 0 ; i < nst->n_vars ; ++i)
235     {
236       const struct mw *mw = &mwv[i];
237       tab_text (table, 0, column_headers + i, TAT_TITLE,
238                 var_to_string (nst->vars[i]));
239
240       tab_double (table, 1, column_headers + i, 0,
241                   mw->n[0], NULL, RC_OTHER);
242
243       tab_double (table, 2, column_headers + i, 0,
244                   mw->n[1], NULL, RC_OTHER);
245
246       tab_double (table, 3, column_headers + i, 0,
247                   mw->n[1] + mw->n[0], NULL, RC_OTHER);
248
249       /* Mean Ranks */
250       tab_double (table, 4, column_headers + i, 0,
251                   mw->rank_sum[0] / mw->n[0], NULL, RC_OTHER);
252
253       tab_double (table, 5, column_headers + i, 0,
254                   mw->rank_sum[1] / mw->n[1], NULL, RC_OTHER);
255
256       /* Sum of Ranks */
257       tab_double (table, 6, column_headers + i, 0,
258                   mw->rank_sum[0], NULL, RC_OTHER);
259
260       tab_double (table, 7, column_headers + i, 0,
261                   mw->rank_sum[1], NULL, RC_OTHER);
262     }
263
264   tab_submit (table);
265 }
266
267 static void
268 show_statistics_box (const struct n_sample_test *nst, const struct mw *mwv, bool exact)
269 {
270   int i;
271   const int row_headers = 1;
272   const int column_headers = 1;
273   struct tab_table *table =
274     tab_create (row_headers + (exact ? 6 : 4), column_headers + nst->n_vars);
275
276   tab_headers (table, row_headers, 0, column_headers, 0);
277
278   tab_title (table, _("Test Statistics"));
279
280   /* Vertical lines inside the box */
281   tab_box (table, 1, 0, -1, TAL_1,
282            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
283
284   /* Box around the table */
285   tab_box (table, TAL_2, TAL_2, -1, -1,
286            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
287
288   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
289   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
290
291   tab_text (table, 1, 0, TAT_TITLE | TAB_CENTER, _("Mann-Whitney U"));
292   tab_text (table, 2, 0, TAT_TITLE | TAB_CENTER, _("Wilcoxon W"));
293   tab_text (table, 3, 0, TAT_TITLE | TAB_CENTER, _("Z"));
294   tab_text (table, 4, 0, TAT_TITLE | TAB_CENTER, _("Asymp. Sig. (2-tailed)"));
295
296   if (exact) 
297     {
298       tab_text (table, 5, 0, TAT_TITLE | TAB_CENTER, _("Exact Sig. (2-tailed)"));
299       tab_text (table, 6, 0, TAT_TITLE | TAB_CENTER, _("Point Probability"));
300     }
301
302   for (i = 0 ; i < nst->n_vars ; ++i)
303     {
304       const struct mw *mw = &mwv[i];
305
306       tab_text (table, 0, column_headers + i, TAT_TITLE,
307                 var_to_string (nst->vars[i]));
308
309       tab_double (table, 1, column_headers + i, 0,
310                   mw->u, NULL, RC_OTHER);
311
312       tab_double (table, 2, column_headers + i, 0,
313                   mw->w, NULL, RC_OTHER);
314
315       tab_double (table, 3, column_headers + i, 0,
316                   mw->z, NULL, RC_OTHER);
317
318       tab_double (table, 4, column_headers + i, 0,
319                   2.0 * gsl_cdf_ugaussian_P (mw->z), NULL, RC_PVALUE);
320     }
321
322   tab_submit (table);
323 }