1 /* Pspp - a program for statistical analysis.
2 Copyright (C) 2008, 2009, 2010, 2011 Free Software Foundation, Inc.
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.
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.
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/>. */
21 #include "language/commands/wilcoxon.h"
23 #include <gsl/gsl_cdf.h>
26 #include "data/casereader.h"
27 #include "data/casewriter.h"
28 #include "data/dataset.h"
29 #include "data/dictionary.h"
30 #include "data/format.h"
31 #include "data/subcase.h"
32 #include "data/variable.h"
33 #include "libpspp/assertion.h"
34 #include "libpspp/message.h"
35 #include "libpspp/misc.h"
36 #include "math/sort.h"
37 #include "math/wilcoxon-sig.h"
38 #include "output/pivot-table.h"
40 #include "gl/minmax.h"
41 #include "gl/xalloc.h"
44 #define N_(msgid) msgid
45 #define _(msgid) gettext (msgid)
48 append_difference (const struct ccase *c, casenumber n UNUSED, void *aux)
50 const variable_pair *vp = aux;
52 return case_num (c, (*vp)[0]) - case_num (c, (*vp)[1]);
55 static void show_ranks_box (const struct wilcoxon_state *,
56 const struct two_sample_test *,
57 const struct dictionary *);
59 static void show_tests_box (const struct wilcoxon_state *,
60 const struct two_sample_test *,
61 bool exact, double timer);
66 distinct_callback (double v UNUSED, casenumber n, double w UNUSED, void *aux)
68 struct wilcoxon_state *ws = aux;
70 ws->tiebreaker += pow3 (n) - n;
76 wilcoxon_execute (const struct dataset *ds,
77 struct casereader *input,
78 enum mv_class exclude,
79 const struct npar_test *test,
85 const struct dictionary *dict = dataset_dict (ds);
86 const struct two_sample_test *t2s = UP_CAST (test, const struct two_sample_test, parent);
88 struct wilcoxon_state *ws = XCALLOC (t2s->n_pairs, struct wilcoxon_state);
89 const struct variable *weight = dict_get_weight (dict);
90 struct caseproto *proto;
93 casereader_create_filter_weight (input, dict, &warn, NULL);
95 proto = caseproto_create ();
96 proto = caseproto_add_width (proto, 0);
97 proto = caseproto_add_width (proto, 0);
99 proto = caseproto_add_width (proto, 0);
101 for (i = 0 ; i < t2s->n_pairs; ++i)
103 struct casereader *r = casereader_clone (input);
104 struct casewriter *writer;
106 struct subcase ordering;
107 variable_pair *vp = &t2s->pairs[i];
109 ws[i].dict = dict_create ("UTF-8");
110 ws[i].sign = dict_create_var (ws[i].dict, "sign", 0);
111 ws[i].absdiff = dict_create_var (ws[i].dict, "absdiff", 0);
112 ws[i].weight = dict_create_var (ws[i].dict, "weight", 0);
114 r = casereader_create_filter_missing (r, *vp, 2,
118 subcase_init_var (&ordering, ws[i].absdiff, SC_ASCEND);
119 writer = sort_create_writer (&ordering, proto);
120 subcase_uninit (&ordering);
122 for (; (c = casereader_read (r)) != NULL; case_unref (c))
124 struct ccase *output = case_create (proto);
125 double d = append_difference (c, 0, vp);
128 *case_num_rw (output, ws[i].sign) = 1.0;
130 *case_num_rw (output, ws[i].sign) = -1.0;
135 w = case_num (c, weight);
137 /* Central point values should be dropped */
143 *case_num_rw (output, ws[i].absdiff) = fabs (d);
146 *case_num_rw (output, ws[i].weight) = case_num (c, weight);
148 casewriter_write (writer, output);
150 casereader_destroy (r);
151 ws[i].reader = casewriter_make_reader (writer);
153 caseproto_unref (proto);
155 for (i = 0 ; i < t2s->n_pairs; ++i)
157 struct casereader *rr ;
159 enum rank_error err = 0;
161 rr = casereader_create_append_rank (ws[i].reader, ws[i].absdiff,
162 weight ? ws[i].weight : NULL, &err,
163 distinct_callback, &ws[i]
166 for (; (c = casereader_read (rr)) != NULL; case_unref (c))
168 double sign = case_num (c, ws[i].sign);
169 double rank = case_num_idx (c, weight ? 3 : 2);
170 double w = weight ? case_num (c, ws[i].weight) : 1.0;
174 ws[i].positives.sum += rank * w;
175 ws[i].positives.n += w;
179 ws[i].negatives.sum += rank * w;
180 ws[i].negatives.n += w;
186 casereader_destroy (rr);
189 casereader_destroy (input);
191 show_ranks_box (ws, t2s, dict);
192 show_tests_box (ws, t2s, exact, timer);
194 for (i = 0 ; i < t2s->n_pairs; ++i)
195 dict_unref (ws[i].dict);
201 put_row (struct pivot_table *table, int var_idx, int sign_idx,
202 double n, double sum)
204 pivot_table_put3 (table, 0, sign_idx, var_idx, pivot_value_new_number (n));
207 pivot_table_put3 (table, 1, sign_idx, var_idx,
208 pivot_value_new_number (sum / n));
209 pivot_table_put3 (table, 2, sign_idx, var_idx,
210 pivot_value_new_number (sum));
215 add_pair_leaf (struct pivot_dimension *dimension, variable_pair *pair)
217 char *label = xasprintf ("%s - %s", var_to_string ((*pair)[0]),
218 var_to_string ((*pair)[1]));
219 return pivot_category_create_leaf (
221 pivot_value_new_user_text_nocopy (label));
225 show_ranks_box (const struct wilcoxon_state *ws,
226 const struct two_sample_test *t2s,
227 const struct dictionary *dict)
229 struct pivot_table *table = pivot_table_create (N_("Ranks"));
230 pivot_table_set_weight_var (table, dict_get_weight (dict));
232 pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"),
233 N_("N"), PIVOT_RC_COUNT,
234 N_("Mean Rank"), PIVOT_RC_OTHER,
235 N_("Sum of Ranks"), PIVOT_RC_OTHER);
237 pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Sign"),
238 N_("Negative Ranks"), N_("Positive Ranks"),
239 N_("Ties"), N_("Total"));
241 struct pivot_dimension *pairs = pivot_dimension_create (
242 table, PIVOT_AXIS_ROW, N_("Pairs"));
244 for (size_t i = 0 ; i < t2s->n_pairs; ++i)
246 variable_pair *vp = &t2s->pairs[i];
247 int pair_idx = add_pair_leaf (pairs, vp);
249 const struct wilcoxon_state *w = &ws[i];
250 put_row (table, pair_idx, 0, w->negatives.n, w->negatives.sum);
251 put_row (table, pair_idx, 1, w->positives.n, w->positives.sum);
252 put_row (table, pair_idx, 2, w->n_zeros, SYSMIS);
253 put_row (table, pair_idx, 3,
254 w->n_zeros + w->positives.n + w->negatives.n, SYSMIS);
257 pivot_table_submit (table);
262 show_tests_box (const struct wilcoxon_state *ws,
263 const struct two_sample_test *t2s,
268 struct pivot_table *table = pivot_table_create (N_("Test Statistics"));
270 struct pivot_dimension *statistics = pivot_dimension_create (
271 table, PIVOT_AXIS_ROW, N_("Statistics"),
272 N_("Z"), PIVOT_RC_OTHER,
273 N_("Asymp. Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE);
275 pivot_category_create_leaves (
277 N_("Exact Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE,
278 N_("Exact Sig. (1-tailed)"), PIVOT_RC_SIGNIFICANCE);
280 struct pivot_dimension *pairs = pivot_dimension_create (
281 table, PIVOT_AXIS_COLUMN, N_("Pairs"));
283 struct pivot_footnote *too_many_pairs = pivot_table_create_footnote (
284 table, pivot_value_new_text (
285 N_("Too many pairs to calculate exact significance")));
287 for (size_t i = 0 ; i < t2s->n_pairs; ++i)
289 variable_pair *vp = &t2s->pairs[i];
290 int pair_idx = add_pair_leaf (pairs, vp);
292 double n = ws[i].positives.n + ws[i].negatives.n;
293 double z = MIN (ws[i].positives.sum, ws[i].negatives.sum);
294 z -= n * (n + 1)/ 4.0;
295 z /= sqrt (n * (n + 1) * (2*n + 1)/24.0 - ws[i].tiebreaker / 48.0);
299 entries[n_entries++] = z;
300 entries[n_entries++] = 2.0 * gsl_cdf_ugaussian_P (z);
302 int footnote_idx = -1;
305 double p = LevelOfSignificanceWXMPSR (ws[i].positives.sum, n);
308 footnote_idx = n_entries;
309 entries[n_entries++] = SYSMIS;
313 entries[n_entries++] = p;
314 entries[n_entries++] = p / 2.0;
318 for (int j = 0; j < n_entries; j++)
320 struct pivot_value *value = pivot_value_new_number (entries[j]);
321 if (j == footnote_idx)
322 pivot_value_add_footnote (value, too_many_pairs);
323 pivot_table_put2 (table, j, pair_idx, value);
327 pivot_table_submit (table);