src/language/stats/runs.c

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