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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2010, 2011, 2014 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/cochran.h"
  20
  21 #include <float.h>
  22 #include <gsl/gsl_cdf.h>
  23 #include <stdbool.h>
  24
  25 #include "data/casereader.h"
  26 #include "data/dataset.h"
  27 #include "data/dictionary.h"
  28 #include "data/format.h"
  29 #include "data/val-type.h"
  30 #include "data/variable.h"
  31 #include "language/stats/npar.h"
  32 #include "libpspp/cast.h"
  33 #include "libpspp/message.h"
  34 #include "libpspp/misc.h"
  35 #include "output/tab.h"
  36
  37 #include "gettext.h"
  38 #define _(msgid) gettext (msgid)
  39
  40 struct cochran
  41 {
  42   double success;
  43   double failure;
  44
  45   double *hits;
  46   double *misses;
  47
  48   const struct dictionary *dict;
  49   double cc;
  50   double df;
  51   double q;
  52 };
  53
  54 static void show_freqs_box (const struct one_sample_test *ost, const struct cochran *ch);
  55 static void show_sig_box (const struct cochran *ch);
  56
  57 void
  58 cochran_execute (const struct dataset *ds,
  59               struct casereader *input,
  60               enum mv_class exclude,
  61               const struct npar_test *test,
  62               bool exact UNUSED, double timer UNUSED)
  63 {
  64   struct one_sample_test *ct = UP_CAST (test, struct one_sample_test, parent);
  65   int v;
  66   struct cochran ch;
  67   const struct dictionary *dict = dataset_dict (ds);
  68   const struct variable *weight = dict_get_weight (dict);
  69
  70   struct ccase *c;
  71   double rowsq = 0;
  72   ch.cc = 0.0;
  73   ch.dict = dict;
  74   ch.success = SYSMIS;
  75   ch.failure = SYSMIS;
  76   ch.hits = xcalloc (ct->n_vars, sizeof *ch.hits);
  77   ch.misses = xcalloc (ct->n_vars, sizeof *ch.misses);
  78
  79   for (; (c = casereader_read (input)); case_unref (c))
  80     {
  81       double case_hits = 0.0;
  82       const double w = weight ? case_data (c, weight)->f: 1.0;
  83       for (v = 0; v < ct->n_vars; ++v)
  84         {
  85           const struct variable *var = ct->vars[v];
  86           const union value *val = case_data (c, var);
  87
  88           if ( var_is_value_missing (var, val, exclude))
  89             continue;
  90
  91           if ( ch.success == SYSMIS)
  92             {
  93               ch.success = val->f;
  94             }
  95           else if (ch.failure == SYSMIS && val->f != ch.success)
  96             {
  97               ch.failure = val->f;
  98             }
  99           if ( ch.success == val->f)
 100             {
 101               ch.hits[v] += w;
 102               case_hits += w;
 103             }
 104           else if ( ch.failure == val->f)
 105             {
 106               ch.misses[v] += w;
 107             }
 108           else
 109             {
 110               msg (MW, _("More than two values encountered.  Cochran Q test will not be run."));
 111               goto finish;
 112             }
 113         }
 114       ch.cc += w;
 115       rowsq += pow2 (case_hits);
 116     }
 117   casereader_destroy (input);
 118
 119   {
 120     double c_l = 0;
 121     double c_l2 = 0;
 122     for (v = 0; v < ct->n_vars; ++v)
 123       {
 124         c_l += ch.hits[v];
 125         c_l2 += pow2 (ch.hits[v]);
 126       }
 127
 128     ch.q = ct->n_vars * c_l2;
 129     ch.q -= pow2 (c_l);
 130     ch.q *= ct->n_vars - 1;
 131
 132     ch.q /= ct->n_vars * c_l - rowsq;
 133
 134     ch.df = ct->n_vars - 1;
 135   }
 136
 137   show_freqs_box (ct, &ch);
 138   show_sig_box (&ch);
 139
 140  finish:
 141
 142   free (ch.hits);
 143   free (ch.misses);
 144 }
 145
 146 static void
 147 show_freqs_box (const struct one_sample_test *ost, const struct cochran *ct)
 148 {
 149   int i;
 150   const struct variable *weight = dict_get_weight (ct->dict);
 151   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 152
 153   const int row_headers = 1;
 154   const int column_headers = 2;
 155   struct tab_table *table =
 156     tab_create (row_headers + 2, column_headers + ost->n_vars);
 157   tab_set_format (table, RC_WEIGHT, wfmt);
 158
 159   tab_headers (table, row_headers, 0, column_headers, 0);
 160
 161   tab_title (table, _("Frequencies"));
 162
 163   /* Vertical lines inside the box */
 164   tab_box (table, 1, 0, -1, TAL_1,
 165            row_headers, 0, tab_nc (table) - 1, tab_nr (table) - 1 );
 166
 167   /* Box around the table */
 168   tab_box (table, TAL_2, TAL_2, -1, -1,
 169            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 170
 171   tab_joint_text (table, 1, 0, 2, 0,
 172                   TAT_TITLE | TAB_CENTER, _("Value"));
 173
 174   tab_text_format (table, 1, 1, 0, _("Success (%.*g)"),
 175                    DBL_DIG + 1, ct->success);
 176   tab_text_format (table, 2, 1, 0, _("Failure (%.*g)"),
 177                    DBL_DIG + 1, ct->failure);
 178
 179   tab_hline (table, TAL_2, 0, tab_nc (table) - 1, column_headers);
 180   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 181
 182   for (i = 0 ; i < ost->n_vars ; ++i)
 183     {
 184       tab_text (table, 0, column_headers + i,
 185                 TAB_LEFT, var_to_string (ost->vars[i]));
 186
 187       tab_double (table, 1, column_headers + i, 0,
 188                   ct->hits[i], NULL, RC_WEIGHT);
 189
 190       tab_double (table, 2, column_headers + i, 0,
 191                   ct->misses[i], NULL, RC_WEIGHT);
 192     }
 193
 194   tab_submit (table);
 195 }
 196
 197
 198
 199 static void
 200 show_sig_box (const struct cochran *ch)
 201 {
 202   const struct variable *weight = dict_get_weight (ch->dict);
 203   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 204
 205   const int row_headers = 1;
 206   const int column_headers = 0;
 207   struct tab_table *table =
 208     tab_create (row_headers + 1, column_headers + 4);
 209
 210
 211   tab_set_format (table, RC_WEIGHT, wfmt);
 212
 213   tab_headers (table, row_headers, 0, column_headers, 0);
 214
 215   tab_title (table, _("Test Statistics"));
 216
 217   tab_text (table,  0, column_headers,
 218             TAT_TITLE | TAB_LEFT , _("N"));
 219
 220   tab_text (table,  0, 1 + column_headers,
 221             TAT_TITLE | TAB_LEFT , _("Cochran's Q"));
 222
 223   tab_text (table,  0, 2 + column_headers,
 224             TAT_TITLE | TAB_LEFT, _("df"));
 225
 226   tab_text (table,  0, 3 + column_headers,
 227             TAT_TITLE | TAB_LEFT, _("Asymp. Sig."));
 228
 229   /* Box around the table */
 230   tab_box (table, TAL_2, TAL_2, -1, -1,
 231            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 232
 233   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 234   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 235
 236   tab_double (table, 1, column_headers,
 237               0, ch->cc, NULL, RC_WEIGHT);
 238
 239   tab_double (table, 1, column_headers + 1,
 240               0, ch->q, NULL, RC_OTHER);
 241
 242   tab_double (table, 1, column_headers + 2,
 243               0, ch->df, NULL, RC_INTEGER);
 244
 245   tab_double (table, 1, column_headers + 3,
 246               0, gsl_cdf_chisq_Q (ch->q, ch->df),
 247               NULL, RC_PVALUE);
 248
 249   tab_submit (table);
 250 }