src/language/stats/binomial.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 2006 Free Software Foundation, Inc.
   3
   4    This program is free software; you can redistribute it and/or
   5    modify it under the terms of the GNU General Public License as
   6    published by the Free Software Foundation; either version 2 of the
   7    License, or (at your option) any later version.
   8
   9    This program is distributed in the hope that it will be useful, but
  10    WITHOUT ANY WARRANTY; without even the implied warranty of
  11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  12    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, write to the Free Software
  16    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  17    02110-1301, USA. */
  18
  19 #include <config.h>
  20 #include <libpspp/compiler.h>
  21 #include <output/table.h>
  22 #include <libpspp/alloc.h>
  23
  24 #include <data/case.h>
  25 #include <data/casefile.h>
  26 #include <data/dictionary.h>
  27 #include <data/procedure.h>
  28 #include <data/variable.h>
  29 #include <data/value.h>
  30 #include <data/value-labels.h>
  31 #include <data/casefilter.h>
  32
  33 #include <libpspp/message.h>
  34 #include <libpspp/assertion.h>
  35
  36 #include "binomial.h"
  37 #include "freq.h"
  38
  39 #include "gettext.h"
  40 #define _(msgid) gettext (msgid)
  41
  42 #include <libpspp/misc.h>
  43
  44 #include <gsl/gsl_cdf.h>
  45 #include <gsl/gsl_randist.h>
  46 #include <gsl-extras/gsl-extras.h>
  47
  48 #include <minmax.h>
  49
  50 #include <libpspp/hash.h>
  51
  52 static double calculate_binomial_internal (double n1, double n2,
  53                                            double p);
  54
  55
  56 static void
  57 swap (double *i1, double *i2)
  58 {
  59   double temp = *i1;
  60   *i1 = *i2;
  61   *i2 = temp;
  62 }
  63
  64 static double
  65 calculate_binomial (double n1, double n2, double p)
  66 {
  67   const double n = n1 + n2;
  68   const bool test_reversed = (n1 / n > p ) ;
  69   if ( test_reversed )
  70     {
  71       p = 1 - p ;
  72       swap (&n1, &n2);
  73     }
  74
  75   return calculate_binomial_internal (n1, n2, p);
  76 }
  77
  78 static double
  79 calculate_binomial_internal (double n1, double n2, double p)
  80 {
  81   /* SPSS Statistical Algorithms has completely different and WRONG
  82      advice here. */
  83
  84   double sig1tailed = gslextras_cdf_binomial_P (n1, n1 + n2, p);
  85
  86   if ( p == 0.5 )
  87     return sig1tailed > 0.5 ? 1.0 :sig1tailed * 2.0;
  88
  89   return sig1tailed ;
  90 }
  91
  92 static void
  93 do_binomial (const struct dictionary *dict,
  94              const struct casefile *cf,
  95              const struct binomial_test *bst,
  96              struct freq *cat1,
  97              struct freq *cat2,
  98              const struct casefilter *filter
  99              )
 100 {
 101   bool warn = true;
 102
 103   const struct one_sample_test *ost = (const struct one_sample_test *) bst;
 104   struct ccase c;
 105   struct casereader *r = casefile_get_reader (cf, NULL);
 106
 107   while (casereader_read(r, &c))
 108     {
 109       int v;
 110       double w =
 111         dict_get_case_weight (dict, &c, &warn);
 112
 113       for (v = 0 ; v < ost->n_vars ; ++v )
 114         {
 115           const struct variable *var = ost->vars[v];
 116           const union value *value = case_data (&c, var);
 117
 118           if ( casefilter_variable_missing (filter, &c, var))
 119             break;
 120
 121           if ( NULL == cat1[v].value )
 122             {
 123               cat1[v].value = value_dup (value, var_get_width (var));
 124               cat1[v].count = w;
 125             }
 126           else if ( 0 == compare_values (cat1[v].value, value,
 127                                          var_get_width (var)))
 128             cat1[v].count += w;
 129           else if ( NULL == cat2[v].value )
 130             {
 131               cat2[v].value = value_dup (value, var_get_width (var));
 132               cat2[v].count = w;
 133             }
 134           else if ( 0 == compare_values (cat2[v].value, value,
 135                                          var_get_width (var)))
 136             cat2[v].count += w;
 137           else if ( bst->category1 == SYSMIS)
 138             msg (ME, _("Variable %s is not dichotomous"), var_get_name (var));
 139         }
 140
 141       case_destroy (&c);
 142     }
 143   casereader_destroy (r);
 144 }
 145
 146
 147
 148 void
 149 binomial_execute (const struct dataset *ds,
 150                   const struct casefile *cf,
 151                   struct casefilter *filter,
 152                   const struct npar_test *test)
 153 {
 154   int v;
 155   const struct binomial_test *bst = (const struct binomial_test *) test;
 156   const struct one_sample_test *ost = (const struct one_sample_test*) test;
 157
 158   struct freq *cat1 = xzalloc (sizeof (*cat1) * ost->n_vars);
 159   struct freq *cat2 = xzalloc (sizeof (*cat1) * ost->n_vars);
 160   struct tab_table *table ;
 161
 162   assert ((bst->category1 == SYSMIS) == (bst->category2 == SYSMIS) );
 163
 164   if ( bst->category1 != SYSMIS )
 165     {
 166       union value v;
 167       v.f = bst->category1;
 168       cat1->value = value_dup (&v, 0);
 169     }
 170
 171   if ( bst->category2 != SYSMIS )
 172     {
 173       union value v;
 174       v.f = bst->category2;
 175       cat2->value = value_dup (&v, 0);
 176     }
 177
 178   do_binomial (dataset_dict(ds), cf, bst, cat1, cat2, filter);
 179
 180   table = tab_create (7, ost->n_vars * 3 + 1, 0);
 181
 182   tab_dim (table, tab_natural_dimensions);
 183
 184   tab_title (table, _("Binomial Test"));
 185
 186   tab_headers (table, 2, 0, 1, 0);
 187
 188   tab_box (table, TAL_1, TAL_1, -1, TAL_1,
 189            0, 0, table->nc - 1, tab_nr(table) - 1 );
 190
 191   for (v = 0 ; v < ost->n_vars; ++v)
 192     {
 193       double n_total, sig;
 194       const struct variable *var = ost->vars[v];
 195       tab_hline (table, TAL_1, 0, tab_nc (table) -1, 1 + v * 3);
 196
 197       /* Titles */
 198       tab_text (table, 0, 1 + v * 3, TAB_LEFT,
 199                 var_to_string (var));
 200
 201       tab_text (table, 1, 1 + v * 3, TAB_LEFT,
 202                 _("Group1"));
 203
 204       tab_text (table, 1, 2 + v * 3, TAB_LEFT,
 205                 _("Group2"));
 206
 207       tab_text (table, 1, 3 + v * 3, TAB_LEFT,
 208                 _("Total"));
 209
 210       /* Test Prop */
 211       tab_float (table, 5, 1 + v * 3, TAB_NONE, bst->p, 8, 3);
 212
 213       /* Category labels */
 214       tab_text (table, 2, 1 + v * 3, TAB_NONE,
 215                 var_get_value_name (var, cat1[v].value));
 216
 217       tab_text (table, 2, 2 + v * 3, TAB_NONE,
 218                 var_get_value_name (var, cat2[v].value));
 219
 220       /* Observed N */
 221       tab_float (table, 3, 1 + v * 3, TAB_NONE,
 222                  cat1[v].count, 8, 0);
 223
 224       tab_float (table, 3, 2 + v * 3, TAB_NONE,
 225                  cat2[v].count, 8, 0);
 226
 227       n_total = cat1[v].count + cat2[v].count;
 228
 229
 230       tab_float (table, 3, 3 + v * 3, TAB_NONE,
 231                  n_total, 8, 0);
 232
 233       /* Observed Proportions */
 234
 235       tab_float (table, 4, 1 + v * 3, TAB_NONE,
 236                  cat1[v].count / n_total, 8, 3);
 237
 238       tab_float (table, 4, 2 + v * 3, TAB_NONE,
 239                  cat2[v].count / n_total, 8, 3);
 240
 241       tab_float (table, 4, 3 + v * 3, TAB_NONE,
 242                  (cat1[v].count + cat2[v].count) / n_total, 8, 2);
 243
 244
 245       /* Significance */
 246       sig = calculate_binomial (cat1[v].count, cat2[v].count,
 247                                        bst->p);
 248
 249       tab_float (table, 6, 1 + v * 3, TAB_NONE,
 250                  sig, 8, 3);
 251     }
 252
 253   tab_text (table,  2, 0,  TAB_CENTER, _("Category"));
 254   tab_text (table,  3, 0,  TAB_CENTER, _("N"));
 255   tab_text (table,  4, 0,  TAB_CENTER, _("Observed Prop."));
 256   tab_text (table,  5, 0,  TAB_CENTER, _("Test Prop."));
 257
 258   tab_text (table,  6, 0,  TAB_CENTER | TAT_PRINTF,
 259             _("Exact Sig. (%d-tailed)"),
 260             bst->p == 0.5 ? 2: 1);
 261
 262   tab_vline (table, TAL_2, 2, 0, tab_nr (table) -1);
 263
 264   free (cat1);
 265   free (cat2);
 266
 267   tab_submit (table);
 268
 269 }