X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Fwilcoxon.c;h=e1fb14f01506035d40bcee4f37a9afd17f4788c9;hb=2031b99a3f1970a0b9a840652f0aff80ec34b433;hp=7552d7ece1f84e37acd189de6aaeebf4ce11ba2a;hpb=c3ac5a8af9c449072c7e872ca70a78c1755ae309;p=pspp diff --git a/src/language/stats/wilcoxon.c b/src/language/stats/wilcoxon.c index 7552d7ece1..e1fb14f015 100644 --- a/src/language/stats/wilcoxon.c +++ b/src/language/stats/wilcoxon.c @@ -1,5 +1,5 @@ /* Pspp - a program for statistical analysis. - Copyright (C) 2008 Free Software Foundation, Inc. + Copyright (C) 2008, 2009, 2010, 2011 Free Software Foundation, Inc. This program is free software: you can redistribute it and/or modify it under the terms of the GNU General Public License as published by @@ -17,36 +17,44 @@ #include -#include "wilcoxon.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -static double timed_wilcoxon_significance (double w, long int n, double timer); +#include "language/stats/wilcoxon.h" +#include +#include + +#include "data/casereader.h" +#include "data/casewriter.h" +#include "data/dataset.h" +#include "data/dictionary.h" +#include "data/format.h" +#include "data/subcase.h" +#include "data/variable.h" +#include "libpspp/assertion.h" +#include "libpspp/message.h" +#include "libpspp/misc.h" +#include "math/sort.h" +#include "math/wilcoxon-sig.h" +#include "output/pivot-table.h" + +#include "gl/minmax.h" +#include "gl/xalloc.h" + +#include "gettext.h" +#define N_(msgid) msgid +#define _(msgid) gettext (msgid) static double append_difference (const struct ccase *c, casenumber n UNUSED, void *aux) { const variable_pair *vp = aux; - return case_data (c, (*vp)[0])->f - case_data (c, (*vp)[1])->f; + return case_num (c, (*vp)[0]) - case_num (c, (*vp)[1]); } static void show_ranks_box (const struct wilcoxon_state *, - const struct two_sample_test *); + const struct two_sample_test *, + const struct dictionary *); static void show_tests_box (const struct wilcoxon_state *, const struct two_sample_test *, @@ -75,96 +83,92 @@ wilcoxon_execute (const struct dataset *ds, int i; bool warn = true; const struct dictionary *dict = dataset_dict (ds); - const struct two_sample_test *t2s = (struct two_sample_test *) test; + const struct two_sample_test *t2s = UP_CAST (test, const struct two_sample_test, parent); - struct wilcoxon_state *ws = xcalloc (sizeof (*ws), t2s->n_pairs); + struct wilcoxon_state *ws = XCALLOC (t2s->n_pairs, struct wilcoxon_state); const struct variable *weight = dict_get_weight (dict); - struct variable *weightx = var_create_internal (WEIGHT_IDX); + struct variable *weightx = dict_create_internal_var (WEIGHT_IDX, 0); + struct caseproto *proto; input = casereader_create_filter_weight (input, dict, &warn, NULL); - for (i = 0 ; i < t2s->n_pairs; ++i ) + proto = caseproto_create (); + proto = caseproto_add_width (proto, 0); + proto = caseproto_add_width (proto, 0); + if (weight != NULL) + proto = caseproto_add_width (proto, 0); + + for (i = 0 ; i < t2s->n_pairs; ++i) { struct casereader *r = casereader_clone (input); struct casewriter *writer; - struct ccase c; - struct case_ordering *ordering = case_ordering_create (); + struct ccase *c; + struct subcase ordering; variable_pair *vp = &t2s->pairs[i]; - const int reader_width = weight ? 3 : 2; - - ws[i].sign = var_create_internal (0); - ws[i].absdiff = var_create_internal (1); - - case_ordering_add_var (ordering, ws[i].absdiff, SRT_ASCEND); - + ws[i].sign = dict_create_internal_var (0, 0); + ws[i].absdiff = dict_create_internal_var (1, 0); r = casereader_create_filter_missing (r, *vp, 2, exclude, NULL, NULL); - writer = sort_create_writer (ordering, reader_width); - while (casereader_read (r, &c)) - { - struct ccase output; - double d = append_difference (&c, 0, vp); + subcase_init_var (&ordering, ws[i].absdiff, SC_ASCEND); + writer = sort_create_writer (&ordering, proto); + subcase_uninit (&ordering); - case_create (&output, reader_width); + for (; (c = casereader_read (r)) != NULL; case_unref (c)) + { + struct ccase *output = case_create (proto); + double d = append_difference (c, 0, vp); if (d > 0) - { - case_data_rw (&output, ws[i].sign)->f = 1.0; - - } + *case_num_rw (output, ws[i].sign) = 1.0; else if (d < 0) - { - case_data_rw (&output, ws[i].sign)->f = -1.0; - } + *case_num_rw (output, ws[i].sign) = -1.0; else { double w = 1.0; if (weight) - w = case_data (&c, weight)->f; + w = case_num (c, weight); /* Central point values should be dropped */ ws[i].n_zeros += w; - case_destroy (&c); - continue; + case_unref (output); + continue; } - case_data_rw (&output, ws[i].absdiff)->f = fabs (d); + *case_num_rw (output, ws[i].absdiff) = fabs (d); if (weight) - case_data_rw (&output, weightx)->f = case_data (&c, weight)->f; + *case_num_rw (output, weightx) = case_num (c, weight); - casewriter_write (writer, &output); - case_destroy (&c); + casewriter_write (writer, output); } casereader_destroy (r); ws[i].reader = casewriter_make_reader (writer); } + caseproto_unref (proto); - for (i = 0 ; i < t2s->n_pairs; ++i ) + for (i = 0 ; i < t2s->n_pairs; ++i) { struct casereader *rr ; - struct ccase c; + struct ccase *c; enum rank_error err = 0; rr = casereader_create_append_rank (ws[i].reader, ws[i].absdiff, weight ? weightx : NULL, &err, distinct_callback, &ws[i] - ); + ); - while (casereader_read (rr, &c)) + for (; (c = casereader_read (rr)) != NULL; case_unref (c)) { - double sign = case_data (&c, ws[i].sign)->f; - double rank = case_data_idx (&c, weight ? 3 : 2)->f; - double w = 1.0; - if (weight) - w = case_data (&c, weightx)->f; + double sign = case_num (c, ws[i].sign); + double rank = case_num_idx (c, weight ? 3 : 2); + double w = weight ? case_num (c, weightx) : 1.0; - if ( sign > 0 ) + if (sign > 0) { ws[i].positives.sum += rank * w; ws[i].positives.n += w; @@ -176,8 +180,6 @@ wilcoxon_execute (const struct dataset *ds, } else NOT_REACHED (); - - case_destroy (&c); } casereader_destroy (rr); @@ -185,100 +187,78 @@ wilcoxon_execute (const struct dataset *ds, casereader_destroy (input); - var_destroy (weightx); + dict_destroy_internal_var (weightx); - show_ranks_box (ws, t2s); + show_ranks_box (ws, t2s, dict); show_tests_box (ws, t2s, exact, timer); - for (i = 0 ; i < t2s->n_pairs; ++i ) + for (i = 0 ; i < t2s->n_pairs; ++i) { - var_destroy (ws[i].sign); - var_destroy (ws[i].absdiff); + dict_destroy_internal_var (ws[i].sign); + dict_destroy_internal_var (ws[i].absdiff); } free (ws); } - - - -#include "gettext.h" -#define _(msgid) gettext (msgid) - static void -show_ranks_box (const struct wilcoxon_state *ws, const struct two_sample_test *t2s) +put_row (struct pivot_table *table, int var_idx, int sign_idx, + double n, double sum) { - size_t i; - struct tab_table *table = tab_create (5, 1 + 4 * t2s->n_pairs, 0); - - tab_dim (table, tab_natural_dimensions); - - tab_title (table, _("Ranks")); - - tab_headers (table, 2, 0, 1, 0); + pivot_table_put3 (table, 0, sign_idx, var_idx, pivot_value_new_number (n)); + if (sum != SYSMIS) + { + pivot_table_put3 (table, 1, sign_idx, var_idx, + pivot_value_new_number (sum / n)); + pivot_table_put3 (table, 2, sign_idx, var_idx, + pivot_value_new_number (sum)); + } +} - /* Vertical lines inside the box */ - tab_box (table, 0, 0, -1, TAL_1, - 1, 0, table->nc - 1, tab_nr (table) - 1 ); +static int +add_pair_leaf (struct pivot_dimension *dimension, variable_pair *pair) +{ + char *label = xasprintf ("%s - %s", var_to_string ((*pair)[0]), + var_to_string ((*pair)[1])); + return pivot_category_create_leaf ( + dimension->root, + pivot_value_new_user_text_nocopy (label)); +} - /* Box around entire table */ - tab_box (table, TAL_2, TAL_2, -1, -1, - 0, 0, table->nc - 1, tab_nr (table) - 1 ); +static void +show_ranks_box (const struct wilcoxon_state *ws, + const struct two_sample_test *t2s, + const struct dictionary *dict) +{ + struct pivot_table *table = pivot_table_create (N_("Ranks")); + pivot_table_set_weight_var (table, dict_get_weight (dict)); + pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Statistics"), + N_("N"), PIVOT_RC_COUNT, + N_("Mean Rank"), PIVOT_RC_OTHER, + N_("Sum of Ranks"), PIVOT_RC_OTHER); - tab_text (table, 2, 0, TAB_CENTER, _("N")); - tab_text (table, 3, 0, TAB_CENTER, _("Mean Rank")); - tab_text (table, 4, 0, TAB_CENTER, _("Sum of Ranks")); + pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Sign"), + N_("Negative Ranks"), N_("Positive Ranks"), + N_("Ties"), N_("Total")); + struct pivot_dimension *pairs = pivot_dimension_create ( + table, PIVOT_AXIS_ROW, N_("Pairs")); - for (i = 0 ; i < t2s->n_pairs; ++i) + for (size_t i = 0 ; i < t2s->n_pairs; ++i) { variable_pair *vp = &t2s->pairs[i]; - - struct string pair_name; - ds_init_cstr (&pair_name, var_to_string ((*vp)[0])); - ds_put_cstr (&pair_name, " - "); - ds_put_cstr (&pair_name, var_to_string ((*vp)[1])); - - tab_text (table, 1, 1 + i * 4, TAB_LEFT, _("Negative Ranks")); - tab_text (table, 1, 2 + i * 4, TAB_LEFT, _("Positive Ranks")); - tab_text (table, 1, 3 + i * 4, TAB_LEFT, _("Ties")); - tab_text (table, 1, 4 + i * 4, TAB_LEFT, _("Total")); - - tab_hline (table, TAL_1, 0, table->nc - 1, 1 + i * 4); - - - tab_text (table, 0, 1 + i * 4, TAB_LEFT, ds_cstr (&pair_name)); - ds_destroy (&pair_name); - - - /* N */ - tab_float (table, 2, 1 + i * 4, TAB_RIGHT, ws[i].negatives.n, 8, 0); - tab_float (table, 2, 2 + i * 4, TAB_RIGHT, ws[i].positives.n, 8, 0); - tab_float (table, 2, 3 + i * 4, TAB_RIGHT, ws[i].n_zeros, 8, 0); - - tab_float (table, 2, 4 + i * 4, TAB_RIGHT, - ws[i].n_zeros + ws[i].positives.n + ws[i].negatives.n, 8, 0); - - /* Sums */ - tab_float (table, 4, 1 + i * 4, TAB_RIGHT, ws[i].negatives.sum, 8, 2); - tab_float (table, 4, 2 + i * 4, TAB_RIGHT, ws[i].positives.sum, 8, 2); - - - /* Means */ - tab_float (table, 3, 1 + i * 4, TAB_RIGHT, - ws[i].negatives.sum / (double) ws[i].negatives.n, 8, 2); - - tab_float (table, 3, 2 + i * 4, TAB_RIGHT, - ws[i].positives.sum / (double) ws[i].positives.n, 8, 2); - + int pair_idx = add_pair_leaf (pairs, vp); + + const struct wilcoxon_state *w = &ws[i]; + put_row (table, pair_idx, 0, w->negatives.n, w->negatives.sum); + put_row (table, pair_idx, 1, w->positives.n, w->positives.sum); + put_row (table, pair_idx, 2, w->n_zeros, SYSMIS); + put_row (table, pair_idx, 3, + w->n_zeros + w->positives.n + w->negatives.n, SYSMIS); } - tab_hline (table, TAL_2, 0, table->nc - 1, 1); - tab_vline (table, TAL_2, 2, 0, table->nr - 1); - - - tab_submit (table); + pivot_table_submit (table); } @@ -286,133 +266,67 @@ static void show_tests_box (const struct wilcoxon_state *ws, const struct two_sample_test *t2s, bool exact, - double timer + double timer UNUSED ) { - size_t i; - struct tab_table *table = tab_create (1 + t2s->n_pairs, exact ? 5 : 3, 0); - - tab_dim (table, tab_natural_dimensions); - - tab_title (table, _("Test Statistics")); - - tab_headers (table, 1, 0, 1, 0); - - /* Vertical lines inside the box */ - tab_box (table, 0, 0, -1, TAL_1, - 0, 0, table->nc - 1, tab_nr (table) - 1 ); - - /* Box around entire table */ - tab_box (table, TAL_2, TAL_2, -1, -1, - 0, 0, table->nc - 1, tab_nr (table) - 1 ); - - - tab_text (table, 0, 1, TAB_LEFT, _("Z")); - tab_text (table, 0, 2, TAB_LEFT, _("Asymp. Sig (2-tailed)")); - - if ( exact ) - { - tab_text (table, 0, 3, TAB_LEFT, _("Exact Sig (2-tailed)")); - tab_text (table, 0, 4, TAB_LEFT, _("Exact Sig (1-tailed)")); - -#if 0 - tab_text (table, 0, 5, TAB_LEFT, _("Point Probability")); -#endif - } - - for (i = 0 ; i < t2s->n_pairs; ++i) + struct pivot_table *table = pivot_table_create (N_("Test Statistics")); + + struct pivot_dimension *statistics = pivot_dimension_create ( + table, PIVOT_AXIS_ROW, N_("Statistics"), + N_("Z"), PIVOT_RC_OTHER, + N_("Asymp. Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE); + if (exact) + pivot_category_create_leaves ( + statistics->root, + N_("Exact Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE, + N_("Exact Sig. (1-tailed)"), PIVOT_RC_SIGNIFICANCE); + + struct pivot_dimension *pairs = pivot_dimension_create ( + table, PIVOT_AXIS_COLUMN, N_("Pairs")); + + struct pivot_footnote *too_many_pairs = pivot_table_create_footnote ( + table, pivot_value_new_text ( + N_("Too many pairs to calculate exact significance"))); + + for (size_t i = 0 ; i < t2s->n_pairs; ++i) { - double z; - double n = ws[i].positives.n + ws[i].negatives.n; variable_pair *vp = &t2s->pairs[i]; + int pair_idx = add_pair_leaf (pairs, vp); - struct string pair_name; - ds_init_cstr (&pair_name, var_to_string ((*vp)[0])); - ds_put_cstr (&pair_name, " - "); - ds_put_cstr (&pair_name, var_to_string ((*vp)[1])); - - - tab_text (table, 1 + i, 0, TAB_CENTER, ds_cstr (&pair_name)); - ds_destroy (&pair_name); - - z = MIN (ws[i].positives.sum, ws[i].negatives.sum); + double n = ws[i].positives.n + ws[i].negatives.n; + double z = MIN (ws[i].positives.sum, ws[i].negatives.sum); z -= n * (n + 1)/ 4.0; - z /= sqrt (n * (n + 1) * (2*n + 1)/24.0 - ws[i].tiebreaker / 48.0); - tab_float (table, 1 + i, 1, TAB_RIGHT, z, 8, 3); - - tab_float (table, 1 + i, 2, TAB_RIGHT, - 2.0 * gsl_cdf_ugaussian_P (z), - 8, 3); + double entries[4]; + int n_entries = 0; + entries[n_entries++] = z; + entries[n_entries++] = 2.0 * gsl_cdf_ugaussian_P (z); + int footnote_idx = -1; if (exact) { - double p = - timed_wilcoxon_significance (ws[i].positives.sum, - n, - timer ); - - if ( p == SYSMIS) + double p = LevelOfSignificanceWXMPSR (ws[i].positives.sum, n); + if (p < 0) { - msg (MW, _("Exact significance was not calculated after %.2f minutes. Skipping test."), timer); + footnote_idx = n_entries; + entries[n_entries++] = SYSMIS; } else - { - tab_float (table, 1 + i, 3, TAB_RIGHT, p, 8, 3); - tab_float (table, 1 + i, 4, TAB_RIGHT, p / 2.0, 8, 3); - } - } + { + entries[n_entries++] = p; + entries[n_entries++] = p / 2.0; + } + } + + for (int j = 0; j < n_entries; j++) + { + struct pivot_value *value = pivot_value_new_number (entries[j]); + if (j == footnote_idx) + pivot_value_add_footnote (value, too_many_pairs); + pivot_table_put2 (table, j, pair_idx, value); + } } - tab_hline (table, TAL_2, 0, table->nc - 1, 1); - tab_vline (table, TAL_2, 1, 0, table->nr - 1); - - - tab_submit (table); -} - - - -#include - -static sigjmp_buf env; - -static void -give_up_callback (int signal UNUSED) -{ - siglongjmp (env, 1); -} - -static double -timed_wilcoxon_significance (double w, long int n, double timer) -{ - double p = SYSMIS; - - sigset_t set; - - struct sigaction timeout_action; - struct sigaction old_action; - - if (timer <= 0 ) - return LevelOfSignificanceWXMPSR (w, n); - - sigemptyset (&set); - - timeout_action.sa_mask = set; - timeout_action.sa_flags = 0; - - timeout_action.sa_handler = give_up_callback; - - if ( 0 == sigsetjmp (env, 1)) - { - sigaction (SIGALRM, &timeout_action, &old_action); - alarm (timer * 60.0); - - p = LevelOfSignificanceWXMPSR (w, n); - } - - sigaction (SIGALRM, &old_action, NULL); - - return p; + pivot_table_submit (table); }