src/language/stats/runs.c

   1 /* PSPP - a program for statistical analysis. -*-c-*-
   2    Copyright (C) 2010 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
  18 #include <config.h>
  19
  20 #include "runs.h"
  21
  22 #include <gsl/gsl_cdf.h>
  23 #include <math.h>
  24
  25 #include <data/format.h>
  26
  27 #include <libpspp/misc.h>
  28 #include <libpspp/message.h>
  29 #include <data/procedure.h>
  30 #include <data/casereader.h>
  31 #include <data/casewriter.h>
  32 #include <data/casegrouper.h>
  33 #include <data/dictionary.h>
  34 #include <data/subcase.h>
  35 #include <data/variable.h>
  36 #include <math/percentiles.h>
  37 #include <math/sort.h>
  38
  39
  40 #include "gettext.h"
  41 #define _(msgid) gettext (msgid)
  42
  43
  44 struct run_state
  45 {
  46   /* The value used to dichotimise the data */
  47   double cutpoint;
  48
  49   /* The number of cases not less than cutpoint */
  50   double np;
  51
  52   /* The number of cases less than cutpoint */
  53   double nn;
  54
  55   /* The sum of np and nn */
  56   double n;
  57
  58   /* The number of runs */
  59   long runs;
  60
  61   /* The sign of the last case seen */
  62   short last_sign;
  63 };
  64
  65
  66
  67 /* Return the Z statistic representing the assympototic
  68    distribution of the the number of runs */
  69 static double
  70 runs_statistic (const struct run_state *rs)
  71 {
  72   double z;
  73   double sigma;
  74   double mu  = 2 * rs->np * rs->nn;
  75   mu /= rs->np + rs->nn;
  76   mu += 1.0;
  77
  78   z = rs->runs - mu;
  79
  80   if ( rs->n < 50)
  81     {
  82       if (z <= -0.5)
  83         z += 0.5;
  84       else if (z >= 0.5)
  85         z -= 0.5;
  86       else
  87         return 0;
  88     }
  89
  90   sigma = 2 * rs->np * rs->nn;
  91   sigma *= 2 * rs->np * rs->nn - rs->nn - rs->np;
  92   sigma /= pow2 (rs->np + rs->nn);
  93   sigma /= rs->np + rs->nn - 1.0;
  94   sigma = sqrt (sigma);
  95
  96   z /= sigma;
  97
  98   return z;
  99 }
 100
 101 static void show_runs_result (const struct runs_test *, const struct run_state *, const struct dictionary *);
 102
 103 void
 104 runs_execute (const struct dataset *ds,
 105               struct casereader *input,
 106               enum mv_class exclude,
 107               const struct npar_test *test,
 108               bool exact UNUSED,
 109               double timer UNUSED)
 110 {
 111   int v;
 112   struct ccase *c;
 113   const struct dictionary *dict = dataset_dict (ds);
 114   const struct variable *weight = dict_get_weight (dict);
 115
 116   struct one_sample_test *otp = UP_CAST (test, struct one_sample_test, parent);
 117   struct runs_test *rt = UP_CAST (otp, struct runs_test, parent);
 118   struct run_state *rs = xcalloc (otp->n_vars, sizeof (*rs));
 119
 120   switch  ( rt->cp_mode)
 121     {
 122     case CP_MODE:
 123       {
 124         for (v = 0; v < otp->n_vars; ++v)
 125           {
 126             bool multimodal = false;
 127             struct run_state *run = &rs[v];
 128             double last_cc;
 129             struct casereader *group = NULL;
 130             struct casegrouper *grouper;
 131             struct casereader *reader = casereader_clone (input);
 132             const struct variable *var = otp->vars[v];
 133
 134             reader = sort_execute_1var (reader, var);
 135
 136             grouper = casegrouper_create_vars (reader, &var, 1);
 137             last_cc = SYSMIS;
 138             while (casegrouper_get_next_group (grouper, &group))
 139               {
 140                 double x = SYSMIS;
 141                 double cc = 0.0;
 142                 struct ccase *c;
 143                 for (; (c = casereader_read (group)); case_unref (c))
 144                   {
 145                     const double w = weight ? case_data (c, weight)->f: 1.0;
 146                     const union value *val = case_data (c, var);
 147                     if ( var_is_value_missing (var, val, exclude))
 148                       continue;
 149                     x = val->f;
 150                     cc += w;
 151                   }
 152
 153                 if ( cc > last_cc)
 154                   {
 155                     run->cutpoint = x;
 156                   }
 157                 else if ( cc == last_cc)
 158                   {
 159                     multimodal = true;
 160                     if ( x > run->cutpoint)
 161                       run->cutpoint = x;
 162                   }
 163                 last_cc = cc;
 164                 casereader_destroy (group);
 165               }
 166             casegrouper_destroy (grouper);
 167             if (multimodal)
 168               msg (MW, _("Multiple modes exist for varible `%s'.  Using %g as the threshold value."),
 169                    var_get_name (var), run->cutpoint);
 170           }
 171       }
 172       break;
 173     case CP_MEDIAN:
 174       {
 175         for (v = 0; v < otp->n_vars; ++v)
 176           {
 177             double cc = 0.0;
 178             struct ccase *c;
 179             struct run_state *run = &rs[v];
 180             struct casereader *reader = casereader_clone (input);
 181             const struct variable *var = otp->vars[v];
 182             struct casewriter *writer;
 183             struct percentile *median;
 184             struct order_stats *os;
 185             struct subcase sc;
 186             subcase_init_var (&sc, var, SC_ASCEND);
 187             writer = sort_create_writer (&sc, casereader_get_proto (reader));
 188
 189             for (; (c = casereader_read (reader)); )
 190               {
 191                 const union value *val = case_data (c, var);
 192                 const double w = weight ? case_data (c, weight)->f: 1.0;
 193                 if ( var_is_value_missing (var, val, exclude))
 194                   {
 195                     case_unref (c);
 196                     continue;
 197                   }
 198
 199                 cc += w;
 200                 casewriter_write (writer, c);
 201               }
 202             subcase_destroy (&sc);
 203             casereader_destroy (reader);
 204             reader = casewriter_make_reader (writer);
 205
 206             median = percentile_create (0.5, cc);
 207             os = &median->parent;
 208
 209             order_stats_accumulate (&os, 1,
 210                                     reader,
 211                                     weight,
 212                                     var,
 213                                     exclude);
 214
 215             run->cutpoint = percentile_calculate (median, PC_HAVERAGE);
 216             statistic_destroy (&median->parent.parent);
 217           }
 218       }
 219       break;
 220     case CP_MEAN:
 221       {
 222         struct casereader *reader = casereader_clone (input);
 223         for (; (c = casereader_read (reader)); case_unref (c))
 224           {
 225             const double w = weight ? case_data (c, weight)->f: 1.0;
 226             for (v = 0; v < otp->n_vars; ++v)
 227               {
 228                 const struct variable *var = otp->vars[v];
 229                 const union value *val = case_data (c, var);
 230                 const double x = val->f;
 231                 struct run_state *run = &rs[v];
 232
 233                 if ( var_is_value_missing (var, val, exclude))
 234                   continue;
 235
 236                 run->cutpoint += x * w;
 237                 run->n += w;
 238               }
 239           }
 240         casereader_destroy (reader);
 241         for (v = 0; v < otp->n_vars; ++v)
 242           {
 243             struct run_state *run = &rs[v];
 244             run->cutpoint /= run->n;
 245           }
 246       }
 247       break;
 248     case CP_CUSTOM:
 249       {
 250       for (v = 0; v < otp->n_vars; ++v)
 251         {
 252           struct run_state *run = &rs[v];
 253           run->cutpoint = rt->cutpoint;
 254         }
 255       }
 256       break;
 257     }
 258
 259   for (; (c = casereader_read (input)); case_unref (c))
 260     {
 261       const double w = weight ? case_data (c, weight)->f: 1.0;
 262
 263       for (v = 0; v < otp->n_vars; ++v)
 264         {
 265           struct run_state *run = &rs[v];
 266           const struct variable *var = otp->vars[v];
 267           const union value *val = case_data (c, var);
 268           double x = val->f;
 269           double d = x - run->cutpoint;
 270           short sign = 0;
 271
 272           if ( var_is_value_missing (var, val, exclude))
 273             continue;
 274
 275           if (d >= 0)
 276             {
 277               sign = +1;
 278               run->np += w;
 279             }
 280           else
 281             {
 282               sign = -1;
 283               run->nn += w;
 284             }
 285
 286           if (sign != run->last_sign)
 287             run->runs++;
 288
 289           run->last_sign = sign;
 290         }
 291     }
 292   casereader_destroy (input);
 293
 294   for (v = 0; v < otp->n_vars; ++v)
 295     {
 296       struct run_state *run = &rs[v];
 297       run->n = run->np + run->nn;
 298     }
 299
 300   show_runs_result (rt, rs, dict);
 301
 302   free (rs);
 303 }
 304
 305 \f
 306 #include <output/tab.h>
 307
 308 static void
 309 show_runs_result (const struct runs_test *rt, const struct run_state *rs, const struct dictionary *dict)
 310 {
 311   const struct variable *weight = dict_get_weight (dict);
 312   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 313
 314   const struct one_sample_test *otp = &rt->parent;
 315
 316   int i;
 317   const int row_headers = 1;
 318   const int column_headers = 1;
 319   struct tab_table *table =
 320     tab_create (row_headers + otp->n_vars, column_headers + 7);
 321
 322   tab_headers (table, row_headers, 0, column_headers, 0);
 323
 324   tab_title (table, _("Runs Test"));
 325
 326   /* Box around the table and vertical lines inside*/
 327   tab_box (table, TAL_2, TAL_2, -1, TAL_1,
 328            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 329
 330   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 331   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 332
 333   for (i = 0 ; i < otp->n_vars; ++i)
 334     {
 335       const struct run_state *run = &rs[i];
 336
 337       double z = runs_statistic (run);
 338
 339       tab_text (table,  row_headers + i, 0,
 340                 TAT_TITLE | TAB_CENTER ,
 341                 var_to_string (otp->vars[i]));
 342
 343       tab_double (table, row_headers +i, 1, 0,
 344                   run->cutpoint, 0);
 345
 346       tab_double (table, row_headers +i, 2, 0,
 347                   run->nn, wfmt);
 348
 349       tab_double (table, row_headers +i, 3, 0,
 350                   run->np, wfmt);
 351
 352       tab_double (table, row_headers +i, 4, 0,
 353                   run->n, wfmt);
 354
 355       tab_double (table, row_headers +i, 5, 0,
 356                   run->runs, &F_8_0);
 357
 358       tab_double (table, row_headers +i, 6, 0,
 359                   z, 0);
 360
 361       tab_double (table, row_headers +i, 7, 0,
 362                   2.0 * gsl_cdf_ugaussian_P (z), 0);
 363     }
 364
 365   switch  ( rt->cp_mode)
 366     {
 367     case CP_CUSTOM:
 368       tab_text (table,  0, column_headers ,
 369                 TAT_TITLE | TAB_LEFT , _("Test Value"));
 370       break;
 371     case CP_MODE:
 372       tab_text (table,  0, column_headers ,
 373                 TAT_TITLE | TAB_LEFT , _("Test Value (mode)"));
 374       break;
 375     case CP_MEAN:
 376       tab_text (table,  0, column_headers ,
 377                 TAT_TITLE | TAB_LEFT , _("Test Value (mean)"));
 378       break;
 379     case CP_MEDIAN:
 380       tab_text (table,  0, column_headers ,
 381                 TAT_TITLE | TAB_LEFT , _("Test Value (median)"));
 382       break;
 383     }
 384
 385   tab_text (table,  0, column_headers + 1,
 386             TAT_TITLE | TAB_LEFT , _("Cases < Test Value"));
 387
 388   tab_text (table,  0, column_headers + 2,
 389             TAT_TITLE | TAB_LEFT , _("Cases >= Test Value"));
 390
 391   tab_text (table,  0, column_headers + 3,
 392             TAT_TITLE | TAB_LEFT , _("Total Cases"));
 393
 394   tab_text (table,  0, column_headers + 4,
 395             TAT_TITLE | TAB_LEFT , _("Number of Runs"));
 396
 397   tab_text (table,  0, column_headers + 5,
 398             TAT_TITLE | TAB_LEFT , _("Z"));
 399
 400   tab_text (table,  0, column_headers + 6,
 401             TAT_TITLE | TAB_LEFT , _("Asymp. Sig. (2-tailed)"));
 402
 403   tab_submit (table);
 404 }
 405
 406