Move all command implementations into a single 'commands' directory.
[pspp] / src / language / commands / friedman.c
1 /* PSPP - a program for statistical analysis. -*-c-*-
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
18 #include <config.h>
19
20 #include "language/commands/friedman.h"
21
22 #include <gsl/gsl_cdf.h>
23 #include <math.h>
24
25 #include "data/casereader.h"
26 #include "data/dataset.h"
27 #include "data/dictionary.h"
28 #include "data/format.h"
29 #include "data/variable.h"
30 #include "libpspp/message.h"
31 #include "libpspp/misc.h"
32 #include "output/pivot-table.h"
33
34 #include "gettext.h"
35 #define N_(msgid) msgid
36 #define _(msgid) gettext (msgid)
37
38
39 struct friedman
40 {
41   double *rank_sum;
42   double cc;
43   double chi_sq;
44   double w;
45   const struct dictionary *dict;
46 };
47
48 static void show_ranks_box (const struct one_sample_test *ost,
49                             const struct friedman *fr);
50
51 static void show_sig_box (const struct one_sample_test *ost,
52                           const struct friedman *fr);
53
54 struct datum
55 {
56   long posn;
57   double x;
58 };
59
60 static int
61 cmp_x (const void *a_, const void *b_)
62 {
63   const struct datum *a = a_;
64   const struct datum *b = b_;
65
66   if (a->x < b->x)
67     return -1;
68
69   return (a->x > b->x);
70 }
71
72 static int
73 cmp_posn (const void *a_, const void *b_)
74 {
75   const struct datum *a = a_;
76   const struct datum *b = b_;
77
78   if (a->posn < b->posn)
79     return -1;
80
81   return (a->posn > b->posn);
82 }
83
84 void
85 friedman_execute (const struct dataset *ds,
86                   struct casereader *input,
87                   enum mv_class exclude,
88                   const struct npar_test *test,
89                   bool exact UNUSED,
90                   double timer UNUSED)
91 {
92   double numerator = 0.0;
93   double denominator = 0.0;
94   int v;
95   struct ccase *c;
96   const struct dictionary *dict = dataset_dict (ds);
97   const struct variable *weight = dict_get_weight (dict);
98
99   struct one_sample_test *ost = UP_CAST (test, struct one_sample_test, parent);
100   struct friedman_test *ft = UP_CAST (ost, struct friedman_test, parent);
101   bool warn = true;
102
103   double sigma_t = 0.0;
104   struct datum *row = XCALLOC (ost->n_vars,  struct datum);
105   double rsq;
106   struct friedman fr;
107   fr.rank_sum = xcalloc (ost->n_vars, sizeof *fr.rank_sum);
108   fr.cc = 0.0;
109   fr.dict = dict;
110   for (v = 0; v < ost->n_vars; ++v)
111     {
112       row[v].posn = v;
113       fr.rank_sum[v] = 0.0;
114     }
115
116   input = casereader_create_filter_weight (input, dict, &warn, NULL);
117   input = casereader_create_filter_missing (input,
118                                             ost->vars, ost->n_vars,
119                                             exclude, 0, 0);
120
121   for (; (c = casereader_read (input)); case_unref (c))
122     {
123       double prev_x = SYSMIS;
124       int run_length = 0;
125
126       const double w = weight ? case_num (c, weight) : 1.0;
127
128       fr.cc += w;
129
130       for (v = 0; v < ost->n_vars; ++v)
131         {
132           const struct variable *var = ost->vars[v];
133           const union value *val = case_data (c, var);
134           row[v].x = val->f;
135         }
136
137       qsort (row, ost->n_vars, sizeof *row, cmp_x);
138       for (v = 0; v < ost->n_vars; ++v)
139         {
140           double x = row[v].x;
141           /* Replace value by the Rank */
142           if (prev_x == x)
143             {
144               /* Deal with ties */
145               int i;
146               run_length++;
147               for (i = v - run_length; i < v; ++i)
148                 {
149                   row[i].x *= run_length ;
150                   row[i].x += v + 1;
151                   row[i].x /= run_length + 1;
152                 }
153               row[v].x = row[v-1].x;
154             }
155           else
156             {
157               row[v].x = v + 1;
158               if (run_length > 0)
159                 {
160                   double t = run_length + 1;
161                   sigma_t += w * (pow3 (t) - t);
162                 }
163               run_length = 0;
164             }
165           prev_x = x;
166         }
167       if (run_length > 0)
168         {
169           double t = run_length + 1;
170           sigma_t += w * (pow3 (t) - t);
171         }
172
173       qsort (row, ost->n_vars, sizeof *row, cmp_posn);
174
175       for (v = 0; v < ost->n_vars; ++v)
176         fr.rank_sum[v] += row[v].x * w;
177     }
178   casereader_destroy (input);
179   free (row);
180
181
182   for (v = 0; v < ost->n_vars; ++v)
183     {
184       numerator += pow2 (fr.rank_sum[v]);
185     }
186
187   rsq = numerator;
188
189   numerator *= 12.0 / (fr.cc * ost->n_vars * (ost->n_vars + 1));
190   numerator -= 3 * fr.cc * (ost->n_vars + 1);
191
192   denominator = 1 - sigma_t / (fr.cc * ost->n_vars * (pow2 (ost->n_vars) - 1));
193
194   fr.chi_sq = numerator / denominator;
195
196   if (ft->kendalls_w)
197     {
198       fr.w = 12 * rsq ;
199       fr.w -= 3 * pow2 (fr.cc) *
200         ost->n_vars * pow2 (ost->n_vars + 1);
201
202       fr.w /= pow2 (fr.cc) * (pow3 (ost->n_vars) - ost->n_vars)
203         - fr.cc * sigma_t;
204     }
205   else
206     fr.w = SYSMIS;
207
208   show_ranks_box (ost, &fr);
209   show_sig_box (ost, &fr);
210
211   free (fr.rank_sum);
212 }
213
214 \f
215
216
217 static void
218 show_ranks_box (const struct one_sample_test *ost, const struct friedman *fr)
219 {
220   struct pivot_table *table = pivot_table_create (N_("Ranks"));
221
222   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Mean Rank"),
223                           N_("Mean Rank"), PIVOT_RC_OTHER);
224
225   struct pivot_dimension *variables = pivot_dimension_create (
226     table, PIVOT_AXIS_ROW, N_("Variable"));
227
228   for (size_t i = 0 ; i < ost->n_vars ; ++i)
229     {
230       int row = pivot_category_create_leaf (
231         variables->root, pivot_value_new_variable (ost->vars[i]));
232
233       pivot_table_put2 (table, 0, row,
234                         pivot_value_new_number (fr->rank_sum[i] / fr->cc));
235     }
236
237   pivot_table_submit (table);
238 }
239
240
241 static void
242 show_sig_box (const struct one_sample_test *ost, const struct friedman *fr)
243 {
244   const struct friedman_test *ft = UP_CAST (ost, const struct friedman_test, parent);
245
246   struct pivot_table *table = pivot_table_create (N_("Test Statistics"));
247   pivot_table_set_weight_var (table, dict_get_weight (fr->dict));
248
249   struct pivot_dimension *statistics = pivot_dimension_create (
250     table, PIVOT_AXIS_ROW, N_("Statistics"),
251     N_("N"), PIVOT_RC_COUNT);
252   if (ft->kendalls_w)
253     pivot_category_create_leaves (statistics->root, N_("Kendall's W"),
254                                   PIVOT_RC_OTHER);
255   pivot_category_create_leaves (statistics->root,
256                                 N_("Chi-Square"), PIVOT_RC_OTHER,
257                                 N_("df"), PIVOT_RC_INTEGER,
258                                 N_("Asymp. Sig."), PIVOT_RC_SIGNIFICANCE);
259
260   double entries[5];
261   int n = 0;
262
263   entries[n++] = fr->cc;
264   if (ft->kendalls_w)
265     entries[n++] = fr->w;
266   entries[n++] = fr->chi_sq;
267   entries[n++] = ost->n_vars - 1;
268   entries[n++] = gsl_cdf_chisq_Q (fr->chi_sq, ost->n_vars - 1);
269   assert (n <= sizeof entries / sizeof *entries);
270
271   for (size_t i = 0; i < n; i++)
272     pivot_table_put1 (table, i, pivot_value_new_number (entries[i]));
273
274   pivot_table_submit (table);
275 }