X-Git-Url: https://pintos-os.org/cgi-bin/gitweb.cgi?a=blobdiff_plain;f=src%2Flanguage%2Fstats%2Fwilcoxon.c;h=e1fb14f01506035d40bcee4f37a9afd17f4788c9;hb=591bdbcf101c4d28cee791cc2444c414e460f7be;hp=e44233a7ab7904fb437a1988f6505cb3af4559cf;hpb=8953baa61127d6d3b91f763663ea647bf3e4e793;p=pspp diff --git a/src/language/stats/wilcoxon.c b/src/language/stats/wilcoxon.c index e44233a7ab..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, 2009 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,33 +17,44 @@ #include -#include "wilcoxon.h" -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include -#include + +#include "language/stats/wilcoxon.h" + #include -#include -#include -#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 *, @@ -72,16 +83,23 @@ 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; @@ -89,56 +107,51 @@ wilcoxon_execute (const struct dataset *ds, 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); + 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); subcase_init_var (&ordering, ws[i].absdiff, SC_ASCEND); - writer = sort_create_writer (&ordering, reader_width); - subcase_destroy (&ordering); + writer = sort_create_writer (&ordering, proto); + subcase_uninit (&ordering); for (; (c = casereader_read (r)) != NULL; case_unref (c)) { - struct ccase *output = case_create (reader_width); + 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_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); } 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; @@ -147,17 +160,15 @@ wilcoxon_execute (const struct dataset *ds, rr = casereader_create_append_rank (ws[i].reader, ws[i].absdiff, weight ? weightx : NULL, &err, distinct_callback, &ws[i] - ); + ); 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,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); } @@ -280,81 +269,64 @@ show_tests_box (const struct wilcoxon_state *ws, 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 ) + 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) { - 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) - { - 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 = LevelOfSignificanceWXMPSR (ws[i].positives.sum, n); if (p < 0) { - msg (MW, ("Too many pairs to calculate exact significance.")); + 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); + pivot_table_submit (table); }