1 /* Pspp - a program for statistical analysis.
2 Copyright (C) 2008, 2009 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 <data/variable.h>
22 #include <data/casereader.h>
23 #include <data/casewriter.h>
24 #include <data/subcase.h>
25 #include <math/sort.h>
26 #include <libpspp/message.h>
28 #include <output/table.h>
29 #include <data/procedure.h>
30 #include <data/dictionary.h>
31 #include <misc/wx-mp-sr.h>
32 #include <gsl/gsl_cdf.h>
35 #include <libpspp/assertion.h>
37 static double timed_wilcoxon_significance (double w, long int n, double timer);
41 append_difference (const struct ccase *c, casenumber n UNUSED, void *aux)
43 const variable_pair *vp = aux;
45 return case_data (c, (*vp)[0])->f - case_data (c, (*vp)[1])->f;
48 static void show_ranks_box (const struct wilcoxon_state *,
49 const struct two_sample_test *);
51 static void show_tests_box (const struct wilcoxon_state *,
52 const struct two_sample_test *,
53 bool exact, double timer);
58 distinct_callback (double v UNUSED, casenumber n, double w UNUSED, void *aux)
60 struct wilcoxon_state *ws = aux;
62 ws->tiebreaker += pow3 (n) - n;
68 wilcoxon_execute (const struct dataset *ds,
69 struct casereader *input,
70 enum mv_class exclude,
71 const struct npar_test *test,
77 const struct dictionary *dict = dataset_dict (ds);
78 const struct two_sample_test *t2s = (struct two_sample_test *) test;
80 struct wilcoxon_state *ws = xcalloc (sizeof (*ws), t2s->n_pairs);
81 const struct variable *weight = dict_get_weight (dict);
82 struct variable *weightx = var_create_internal (WEIGHT_IDX);
85 casereader_create_filter_weight (input, dict, &warn, NULL);
87 for (i = 0 ; i < t2s->n_pairs; ++i )
89 struct casereader *r = casereader_clone (input);
90 struct casewriter *writer;
92 struct subcase ordering;
93 variable_pair *vp = &t2s->pairs[i];
95 const int reader_width = weight ? 3 : 2;
97 ws[i].sign = var_create_internal (0);
98 ws[i].absdiff = var_create_internal (1);
100 r = casereader_create_filter_missing (r, *vp, 2,
104 subcase_init_var (&ordering, ws[i].absdiff, SC_ASCEND);
105 writer = sort_create_writer (&ordering, reader_width);
106 subcase_destroy (&ordering);
108 for (; (c = casereader_read (r)) != NULL; case_unref (c))
110 struct ccase *output = case_create (reader_width);
111 double d = append_difference (c, 0, vp);
115 case_data_rw (output, ws[i].sign)->f = 1.0;
120 case_data_rw (output, ws[i].sign)->f = -1.0;
126 w = case_data (c, weight)->f;
128 /* Central point values should be dropped */
133 case_data_rw (output, ws[i].absdiff)->f = fabs (d);
136 case_data_rw (output, weightx)->f = case_data (c, weight)->f;
138 casewriter_write (writer, output);
140 casereader_destroy (r);
141 ws[i].reader = casewriter_make_reader (writer);
144 for (i = 0 ; i < t2s->n_pairs; ++i )
146 struct casereader *rr ;
148 enum rank_error err = 0;
150 rr = casereader_create_append_rank (ws[i].reader, ws[i].absdiff,
151 weight ? weightx : NULL, &err,
152 distinct_callback, &ws[i]
155 for (; (c = casereader_read (rr)) != NULL; case_unref (c))
157 double sign = case_data (c, ws[i].sign)->f;
158 double rank = case_data_idx (c, weight ? 3 : 2)->f;
161 w = case_data (c, weightx)->f;
165 ws[i].positives.sum += rank * w;
166 ws[i].positives.n += w;
170 ws[i].negatives.sum += rank * w;
171 ws[i].negatives.n += w;
177 casereader_destroy (rr);
180 casereader_destroy (input);
182 var_destroy (weightx);
184 show_ranks_box (ws, t2s);
185 show_tests_box (ws, t2s, exact, timer);
187 for (i = 0 ; i < t2s->n_pairs; ++i )
189 var_destroy (ws[i].sign);
190 var_destroy (ws[i].absdiff);
200 #define _(msgid) gettext (msgid)
203 show_ranks_box (const struct wilcoxon_state *ws, const struct two_sample_test *t2s)
206 struct tab_table *table = tab_create (5, 1 + 4 * t2s->n_pairs, 0);
208 tab_dim (table, tab_natural_dimensions);
210 tab_title (table, _("Ranks"));
212 tab_headers (table, 2, 0, 1, 0);
214 /* Vertical lines inside the box */
215 tab_box (table, 0, 0, -1, TAL_1,
216 1, 0, table->nc - 1, tab_nr (table) - 1 );
218 /* Box around entire table */
219 tab_box (table, TAL_2, TAL_2, -1, -1,
220 0, 0, table->nc - 1, tab_nr (table) - 1 );
223 tab_text (table, 2, 0, TAB_CENTER, _("N"));
224 tab_text (table, 3, 0, TAB_CENTER, _("Mean Rank"));
225 tab_text (table, 4, 0, TAB_CENTER, _("Sum of Ranks"));
228 for (i = 0 ; i < t2s->n_pairs; ++i)
230 variable_pair *vp = &t2s->pairs[i];
232 struct string pair_name;
233 ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
234 ds_put_cstr (&pair_name, " - ");
235 ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
237 tab_text (table, 1, 1 + i * 4, TAB_LEFT, _("Negative Ranks"));
238 tab_text (table, 1, 2 + i * 4, TAB_LEFT, _("Positive Ranks"));
239 tab_text (table, 1, 3 + i * 4, TAB_LEFT, _("Ties"));
240 tab_text (table, 1, 4 + i * 4, TAB_LEFT, _("Total"));
242 tab_hline (table, TAL_1, 0, table->nc - 1, 1 + i * 4);
245 tab_text (table, 0, 1 + i * 4, TAB_LEFT, ds_cstr (&pair_name));
246 ds_destroy (&pair_name);
250 tab_float (table, 2, 1 + i * 4, TAB_RIGHT, ws[i].negatives.n, 8, 0);
251 tab_float (table, 2, 2 + i * 4, TAB_RIGHT, ws[i].positives.n, 8, 0);
252 tab_float (table, 2, 3 + i * 4, TAB_RIGHT, ws[i].n_zeros, 8, 0);
254 tab_float (table, 2, 4 + i * 4, TAB_RIGHT,
255 ws[i].n_zeros + ws[i].positives.n + ws[i].negatives.n, 8, 0);
258 tab_float (table, 4, 1 + i * 4, TAB_RIGHT, ws[i].negatives.sum, 8, 2);
259 tab_float (table, 4, 2 + i * 4, TAB_RIGHT, ws[i].positives.sum, 8, 2);
263 tab_float (table, 3, 1 + i * 4, TAB_RIGHT,
264 ws[i].negatives.sum / (double) ws[i].negatives.n, 8, 2);
266 tab_float (table, 3, 2 + i * 4, TAB_RIGHT,
267 ws[i].positives.sum / (double) ws[i].positives.n, 8, 2);
271 tab_hline (table, TAL_2, 0, table->nc - 1, 1);
272 tab_vline (table, TAL_2, 2, 0, table->nr - 1);
280 show_tests_box (const struct wilcoxon_state *ws,
281 const struct two_sample_test *t2s,
287 struct tab_table *table = tab_create (1 + t2s->n_pairs, exact ? 5 : 3, 0);
289 tab_dim (table, tab_natural_dimensions);
291 tab_title (table, _("Test Statistics"));
293 tab_headers (table, 1, 0, 1, 0);
295 /* Vertical lines inside the box */
296 tab_box (table, 0, 0, -1, TAL_1,
297 0, 0, table->nc - 1, tab_nr (table) - 1 );
299 /* Box around entire table */
300 tab_box (table, TAL_2, TAL_2, -1, -1,
301 0, 0, table->nc - 1, tab_nr (table) - 1 );
304 tab_text (table, 0, 1, TAB_LEFT, _("Z"));
305 tab_text (table, 0, 2, TAB_LEFT, _("Asymp. Sig (2-tailed)"));
309 tab_text (table, 0, 3, TAB_LEFT, _("Exact Sig (2-tailed)"));
310 tab_text (table, 0, 4, TAB_LEFT, _("Exact Sig (1-tailed)"));
313 tab_text (table, 0, 5, TAB_LEFT, _("Point Probability"));
317 for (i = 0 ; i < t2s->n_pairs; ++i)
320 double n = ws[i].positives.n + ws[i].negatives.n;
321 variable_pair *vp = &t2s->pairs[i];
323 struct string pair_name;
324 ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
325 ds_put_cstr (&pair_name, " - ");
326 ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
329 tab_text (table, 1 + i, 0, TAB_CENTER, ds_cstr (&pair_name));
330 ds_destroy (&pair_name);
332 z = MIN (ws[i].positives.sum, ws[i].negatives.sum);
333 z -= n * (n + 1)/ 4.0;
335 z /= sqrt (n * (n + 1) * (2*n + 1)/24.0 - ws[i].tiebreaker / 48.0);
337 tab_float (table, 1 + i, 1, TAB_RIGHT, z, 8, 3);
339 tab_float (table, 1 + i, 2, TAB_RIGHT,
340 2.0 * gsl_cdf_ugaussian_P (z),
346 timed_wilcoxon_significance (ws[i].positives.sum,
352 msg (MW, _("Exact significance was not calculated after %.2f minutes. Skipping test."), timer);
356 tab_float (table, 1 + i, 3, TAB_RIGHT, p, 8, 3);
357 tab_float (table, 1 + i, 4, TAB_RIGHT, p / 2.0, 8, 3);
362 tab_hline (table, TAL_2, 0, table->nc - 1, 1);
363 tab_vline (table, TAL_2, 1, 0, table->nr - 1);
373 static sigjmp_buf env;
376 give_up_callback (int signal UNUSED)
382 timed_wilcoxon_significance (double w, long int n, double timer)
388 struct sigaction timeout_action;
389 struct sigaction old_action;
392 return LevelOfSignificanceWXMPSR (w, n);
396 timeout_action.sa_mask = set;
397 timeout_action.sa_flags = 0;
399 timeout_action.sa_handler = give_up_callback;
401 if ( 0 == sigsetjmp (env, 1))
403 sigaction (SIGALRM, &timeout_action, &old_action);
404 alarm (timer * 60.0);
406 p = LevelOfSignificanceWXMPSR (w, n);
409 sigaction (SIGALRM, &old_action, NULL);