src/language/stats/sign.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2009 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 #include "sign.h"
  19
  20 #include <xalloc.h>
  21
  22 #include <data/variable.h>
  23 #include <libpspp/str.h>
  24 #include <output/table.h>
  25 #include <gsl/gsl_cdf.h>
  26 #include <gsl/gsl_randist.h>
  27 #include "npar.h"
  28 #include <data/procedure.h>
  29 #include <data/missing-values.h>
  30 #include <data/dictionary.h>
  31 #include <data/casereader.h>
  32 #include <data/format.h>
  33
  34 #include "gettext.h"
  35 #define _(msgid) gettext (msgid)
  36
  37 struct sign_test_params
  38 {
  39   double pos;
  40   double ties;
  41   double neg;
  42
  43   double one_tailed_sig;
  44   double point_prob;
  45 };
  46
  47
  48 static void
  49 output_frequency_table (const struct two_sample_test *t2s,
  50                         const struct sign_test_params *param,
  51                         const struct dictionary *dict)
  52 {
  53   int i;
  54   struct tab_table *table = tab_create (3, 1 + 4 * t2s->n_pairs, 0);
  55
  56   const struct variable *wv = dict_get_weight (dict);
  57   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : & F_8_0;
  58
  59   tab_dim (table, tab_natural_dimensions);
  60
  61   tab_title (table, _("Frequencies"));
  62
  63   tab_headers (table, 2, 0, 1, 0);
  64
  65   /* Vertical lines inside the box */
  66   tab_box (table, 0, 0, -1, TAL_1,
  67            1, 0, table->nc - 1, tab_nr (table) - 1 );
  68
  69   /* Box around entire table */
  70   tab_box (table, TAL_2, TAL_2, -1, -1,
  71            0, 0, table->nc - 1, tab_nr (table) - 1 );
  72
  73   tab_text (table,  2, 0,  TAB_CENTER, _("N"));
  74
  75   for (i = 0 ; i < t2s->n_pairs; ++i)
  76     {
  77       variable_pair *vp = &t2s->pairs[i];
  78
  79       struct string pair_name;
  80       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
  81       ds_put_cstr (&pair_name, " - ");
  82       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
  83
  84       tab_text (table, 0, 1 + i * 4, TAB_LEFT, ds_cstr (&pair_name));
  85
  86       ds_destroy (&pair_name);
  87
  88       tab_hline (table, TAL_1, 0, table->nc - 1, 1 + i * 4);
  89
  90       tab_text (table,  1, 1 + i * 4,  TAB_LEFT, _("Negative Differences"));
  91       tab_text (table,  1, 2 + i * 4,  TAB_LEFT, _("Positive Differences"));
  92       tab_text (table,  1, 3 + i * 4,  TAB_LEFT, _("Ties"));
  93       tab_text (table,  1, 4 + i * 4,  TAB_LEFT, _("Total"));
  94
  95       tab_double (table, 2, 1 + i * 4, TAB_RIGHT, param[i].neg, wfmt);
  96       tab_double (table, 2, 2 + i * 4, TAB_RIGHT, param[i].pos, wfmt);
  97       tab_double (table, 2, 3 + i * 4, TAB_RIGHT, param[i].ties, wfmt);
  98       tab_double (table, 2, 4 + i * 4, TAB_RIGHT,
  99                  param[i].ties + param[i].neg + param[i].pos, wfmt);
 100     }
 101
 102   tab_submit (table);
 103 }
 104
 105 static void
 106 output_statistics_table (const struct two_sample_test *t2s,
 107                          const struct sign_test_params *param)
 108 {
 109   int i;
 110   struct tab_table *table = tab_create (1 + t2s->n_pairs, 4, 0);
 111
 112   tab_dim (table, tab_natural_dimensions);
 113
 114   tab_title (table, _("Test Statistics"));
 115
 116   tab_headers (table, 0, 1,  0, 1);
 117
 118   tab_hline (table, TAL_2, 0, table->nc - 1, 1);
 119   tab_vline (table, TAL_2, 1, 0, table->nr - 1);
 120
 121
 122   /* Vertical lines inside the box */
 123   tab_box (table, -1, -1, -1, TAL_1,
 124            0, 0,
 125            table->nc - 1, tab_nr (table) - 1);
 126
 127   /* Box around entire table */
 128   tab_box (table, TAL_2, TAL_2, -1, -1,
 129            0, 0, table->nc - 1,
 130            tab_nr (table) - 1);
 131
 132   tab_text (table,  0, 1, TAT_TITLE | TAB_LEFT,
 133             _("Exact Sig. (2-tailed)"));
 134
 135   tab_text (table,  0, 2, TAT_TITLE | TAB_LEFT,
 136             _("Exact Sig. (1-tailed)"));
 137
 138   tab_text (table,  0, 3, TAT_TITLE | TAB_LEFT,
 139             _("Point Probability"));
 140
 141   for (i = 0 ; i < t2s->n_pairs; ++i)
 142     {
 143       variable_pair *vp = &t2s->pairs[i];
 144
 145       struct string pair_name;
 146       ds_init_cstr (&pair_name, var_to_string ((*vp)[0]));
 147       ds_put_cstr (&pair_name, " - ");
 148       ds_put_cstr (&pair_name, var_to_string ((*vp)[1]));
 149
 150       tab_text (table,  1 + i, 0, TAB_LEFT, ds_cstr (&pair_name));
 151       ds_destroy (&pair_name);
 152
 153       tab_double (table, 1 + i, 1, TAB_RIGHT,
 154                   param[i].one_tailed_sig * 2, NULL);
 155
 156       tab_double (table, 1 + i, 2, TAB_RIGHT, param[i].one_tailed_sig, NULL);
 157       tab_double (table, 1 + i, 3, TAB_RIGHT, param[i].point_prob, NULL);
 158     }
 159
 160   tab_submit (table);
 161 }
 162
 163 void
 164 sign_execute (const struct dataset *ds,
 165                   struct casereader *input,
 166                   enum mv_class exclude,
 167                   const struct npar_test *test,
 168                   bool exact UNUSED,
 169                   double timer UNUSED)
 170 {
 171   int i;
 172   bool warn = true;
 173   const struct dictionary *dict = dataset_dict (ds);
 174   const struct two_sample_test *t2s = (const struct two_sample_test *) test;
 175   struct ccase *c;
 176
 177   struct sign_test_params *stp = xcalloc (sizeof *stp, t2s->n_pairs);
 178
 179   struct casereader *r = input;
 180
 181   for (; (c = casereader_read (r)) != NULL; case_unref (c))
 182     {
 183       const double weight = dict_get_case_weight (dict, c, &warn);
 184
 185       for (i = 0 ; i < t2s->n_pairs; ++i )
 186         {
 187           variable_pair *vp = &t2s->pairs[i];
 188           const union value *value0 = case_data (c, (*vp)[0]);
 189           const union value *value1 = case_data (c, (*vp)[1]);
 190           const double diff = value0->f - value1->f;
 191
 192           if (var_is_value_missing ((*vp)[0], value0, exclude))
 193             continue;
 194
 195           if (var_is_value_missing ((*vp)[1], value1, exclude))
 196             continue;
 197
 198           if ( diff > 0)
 199             stp[i].pos += weight;
 200           else if (diff < 0)
 201             stp[i].neg += weight;
 202           else
 203             stp[i].ties += weight;
 204         }
 205     }
 206
 207   casereader_destroy (r);
 208
 209   for (i = 0 ; i < t2s->n_pairs; ++i )
 210     {
 211       int r = MIN (stp[i].pos, stp[i].neg);
 212       stp[i].one_tailed_sig = gsl_cdf_binomial_P (r,
 213                                                   0.5,
 214                                                   stp[i].pos + stp[i].neg);
 215
 216       stp[i].point_prob = gsl_ran_binomial_pdf (r, 0.5,
 217                                                 stp[i].pos + stp[i].neg);
 218     }
 219
 220   output_frequency_table (t2s, stp, dict);
 221
 222   output_statistics_table (t2s, stp);
 223
 224   free (stp);
 225 }