pintos-os.org Git - pspp/blob - src/language/stats/sign.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2009, 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 #include <config.h>
  18
  19 #include "language/stats/sign.h"
  20
  21 #include <gsl/gsl_cdf.h>
  22 #include <gsl/gsl_randist.h>
  23
  24 #include "data/casereader.h"
  25 #include "data/dataset.h"
  26 #include "data/dictionary.h"
  27 #include "data/format.h"
  28 #include "data/missing-values.h"
  29 #include "data/variable.h"
  30 #include "language/stats/npar.h"
  31 #include "libpspp/str.h"
  32 #include "output/pivot-table.h"
  33
  34 #include "gl/minmax.h"
  35 #include "gl/xalloc.h"
  36
  37 #include "gettext.h"
  38 #define N_(msgid) msgid
  39 #define _(msgid) gettext (msgid)
  40
  41 struct sign_test_params
  42 {
  43   double pos;
  44   double ties;
  45   double neg;
  46
  47   double one_tailed_sig;
  48   double point_prob;
  49 };
  50
  51 static int
  52 add_pair_leaf (struct pivot_dimension *dimension, variable_pair *pair)
  53 {
  54   char *label = xasprintf ("%s - %s", var_to_string ((*pair)[0]),
  55                            var_to_string ((*pair)[1]));
  56   return pivot_category_create_leaf (
  57     dimension->root,
  58     pivot_value_new_user_text_nocopy (label));
  59 }
  60
  61 static void
  62 output_frequency_table (const struct two_sample_test *t2s,
  63                         const struct sign_test_params *param,
  64                         const struct dictionary *dict)
  65 {
  66   struct pivot_table *table = pivot_table_create (N_("Frequencies"));
  67   pivot_table_set_weight_var (table, dict_get_weight (dict));
  68
  69   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("N"),
  70                           N_("N"), PIVOT_RC_COUNT);
  71
  72   pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Differences"),
  73                           N_("Negative Differences"),
  74                           N_("Positive Differences"),
  75                           N_("Ties"), N_("Total"));
  76
  77   struct pivot_dimension *pairs = pivot_dimension_create (
  78     table, PIVOT_AXIS_ROW, N_("Pairs"));
  79
  80   for (size_t i = 0 ; i < t2s->n_pairs; ++i)
  81     {
  82       variable_pair *vp = &t2s->pairs[i];
  83
  84       int pair_idx = add_pair_leaf (pairs, vp);
  85
  86       const struct sign_test_params *p = &param[i];
  87       double values[] = { p->neg, p->pos, p->ties, p->ties + p->neg + p->pos };
  88       for (size_t j = 0; j < sizeof values / sizeof *values; j++)
  89         pivot_table_put3 (table, 0, j, pair_idx,
  90                           pivot_value_new_number (values[j]));
  91     }
  92
  93   pivot_table_submit (table);
  94 }
  95
  96 static void
  97 output_statistics_table (const struct two_sample_test *t2s,
  98                          const struct sign_test_params *param)
  99 {
 100   struct pivot_table *table = pivot_table_create (N_("Test Statistics"));
 101
 102   pivot_dimension_create (table, PIVOT_AXIS_ROW, N_("Statistics"),
 103                           N_("Exact Sig. (2-tailed)"), PIVOT_RC_SIGNIFICANCE,
 104                           N_("Exact Sig. (1-tailed)"), PIVOT_RC_SIGNIFICANCE,
 105                           N_("Point Probability"), PIVOT_RC_SIGNIFICANCE);
 106
 107   struct pivot_dimension *pairs = pivot_dimension_create (
 108     table, PIVOT_AXIS_COLUMN, N_("Pairs"));
 109
 110   for (size_t i = 0 ; i < t2s->n_pairs; ++i)
 111     {
 112       variable_pair *vp = &t2s->pairs[i];
 113       int pair_idx = add_pair_leaf (pairs, vp);
 114
 115       const struct sign_test_params *p = &param[i];
 116       double values[] = { p->one_tailed_sig * 2,
 117                           p->one_tailed_sig,
 118                           p->point_prob };
 119       for (size_t j = 0; j < sizeof values / sizeof *values; j++)
 120         pivot_table_put2 (table, j, pair_idx,
 121                           pivot_value_new_number (values[j]));
 122     }
 123
 124   pivot_table_submit (table);
 125 }
 126
 127 void
 128 sign_execute (const struct dataset *ds,
 129                   struct casereader *input,
 130                   enum mv_class exclude,
 131                   const struct npar_test *test,
 132                   bool exact UNUSED,
 133                   double timer UNUSED)
 134 {
 135   int i;
 136   bool warn = true;
 137   const struct dictionary *dict = dataset_dict (ds);
 138   const struct two_sample_test *t2s = UP_CAST (test, const struct two_sample_test, parent);
 139   struct ccase *c;
 140
 141   struct sign_test_params *stp = XCALLOC (t2s->n_pairs,  struct sign_test_params);
 142
 143   struct casereader *r = input;
 144
 145   for (; (c = casereader_read (r)) != NULL; case_unref (c))
 146     {
 147       const double weight = dict_get_case_weight (dict, c, &warn);
 148
 149       for (i = 0 ; i < t2s->n_pairs; ++i)
 150         {
 151           variable_pair *vp = &t2s->pairs[i];
 152           const union value *value0 = case_data (c, (*vp)[0]);
 153           const union value *value1 = case_data (c, (*vp)[1]);
 154           const double diff = value0->f - value1->f;
 155
 156           if (var_is_value_missing ((*vp)[0], value0) & exclude)
 157             continue;
 158
 159           if (var_is_value_missing ((*vp)[1], value1) & exclude)
 160             continue;
 161
 162           if (diff > 0)
 163             stp[i].pos += weight;
 164           else if (diff < 0)
 165             stp[i].neg += weight;
 166           else
 167             stp[i].ties += weight;
 168         }
 169     }
 170
 171   casereader_destroy (r);
 172
 173   for (i = 0 ; i < t2s->n_pairs; ++i)
 174     {
 175       int r = MIN (stp[i].pos, stp[i].neg);
 176       stp[i].one_tailed_sig = gsl_cdf_binomial_P (r,
 177                                                   0.5,
 178                                                   stp[i].pos + stp[i].neg);
 179
 180       stp[i].point_prob = gsl_ran_binomial_pdf (r, 0.5,
 181                                                 stp[i].pos + stp[i].neg);
 182     }
 183
 184   output_frequency_table (t2s, stp, dict);
 185
 186   output_statistics_table (t2s, stp);
 187
 188   free (stp);
 189 }