pintos-os.org Git - pspp/blob - 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)); case_unref (c))
 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                   continue;
 195
 196                 cc += w;
 197                 casewriter_write (writer, c);
 198               }
 199             subcase_destroy (&sc);
 200             casereader_destroy (reader);
 201             reader = casewriter_make_reader (writer);
 202
 203             median = percentile_create (0.5, cc);
 204             os = &median->parent;
 205
 206             order_stats_accumulate (&os, 1,
 207                                     reader,
 208                                     weight,
 209                                     var,
 210                                     exclude);
 211
 212             run->cutpoint = percentile_calculate (median, PC_HAVERAGE);
 213             statistic_destroy (&median->parent.parent);
 214           }
 215       }
 216       break;
 217     case CP_MEAN:
 218       {
 219         struct casereader *reader = casereader_clone (input);
 220         for (; (c = casereader_read (reader)); case_unref (c))
 221           {
 222             const double w = weight ? case_data (c, weight)->f: 1.0;
 223             for (v = 0; v < otp->n_vars; ++v)
 224               {
 225                 const struct variable *var = otp->vars[v];
 226                 const union value *val = case_data (c, var);
 227                 const double x = val->f;
 228                 struct run_state *run = &rs[v];
 229
 230                 if ( var_is_value_missing (var, val, exclude))
 231                   continue;
 232
 233                 run->cutpoint += x * w;
 234                 run->n += w;
 235               }
 236           }
 237         casereader_destroy (reader);
 238         for (v = 0; v < otp->n_vars; ++v)
 239           {
 240             struct run_state *run = &rs[v];
 241             run->cutpoint /= run->n;
 242           }
 243       }
 244       break;
 245     case CP_CUSTOM:
 246       {
 247       for (v = 0; v < otp->n_vars; ++v)
 248         {
 249           struct run_state *run = &rs[v];
 250           run->cutpoint = rt->cutpoint;
 251         }
 252       }
 253       break;
 254     }
 255
 256   for (; (c = casereader_read (input)); case_unref (c))
 257     {
 258       const double w = weight ? case_data (c, weight)->f: 1.0;
 259
 260       for (v = 0; v < otp->n_vars; ++v)
 261         {
 262           struct run_state *run = &rs[v];
 263           const struct variable *var = otp->vars[v];
 264           const union value *val = case_data (c, var);
 265           double x = val->f;
 266           double d = x - run->cutpoint;
 267           short sign = 0;
 268
 269           if ( var_is_value_missing (var, val, exclude))
 270             continue;
 271
 272           if (d >= 0)
 273             {
 274               sign = +1;
 275               run->np += w;
 276             }
 277           else
 278             {
 279               sign = -1;
 280               run->nn += w;
 281             }
 282
 283           if (sign != run->last_sign)
 284             run->runs++;
 285
 286           run->last_sign = sign;
 287         }
 288     }
 289   casereader_destroy (input);
 290
 291   for (v = 0; v < otp->n_vars; ++v)
 292     {
 293       struct run_state *run = &rs[v];
 294       run->n = run->np + run->nn;
 295     }
 296
 297   show_runs_result (rt, rs, dict);
 298
 299   free (rs);
 300 }
 301
 302 \f
 303 #include <output/tab.h>
 304
 305 static void
 306 show_runs_result (const struct runs_test *rt, const struct run_state *rs, const struct dictionary *dict)
 307 {
 308   const struct variable *weight = dict_get_weight (dict);
 309   const struct fmt_spec *wfmt = weight ? var_get_print_format (weight) : &F_8_0;
 310
 311   const struct one_sample_test *otp = &rt->parent;
 312
 313   int i;
 314   const int row_headers = 1;
 315   const int column_headers = 1;
 316   struct tab_table *table =
 317     tab_create (row_headers + otp->n_vars, column_headers + 7);
 318
 319   tab_headers (table, row_headers, 0, column_headers, 0);
 320
 321   tab_title (table, _("Runs Test"));
 322
 323   /* Box around the table and vertical lines inside*/
 324   tab_box (table, TAL_2, TAL_2, -1, TAL_1,
 325            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 326
 327   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 328   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 329
 330   for (i = 0 ; i < otp->n_vars; ++i)
 331     {
 332       const struct run_state *run = &rs[i];
 333
 334       double z = runs_statistic (run);
 335
 336       tab_text (table,  row_headers + i, 0,
 337                 TAT_TITLE | TAB_CENTER ,
 338                 var_to_string (otp->vars[i]));
 339
 340       tab_double (table, row_headers +i, 1, 0,
 341                   run->cutpoint, 0);
 342
 343       tab_double (table, row_headers +i, 2, 0,
 344                   run->nn, wfmt);
 345
 346       tab_double (table, row_headers +i, 3, 0,
 347                   run->np, wfmt);
 348
 349       tab_double (table, row_headers +i, 4, 0,
 350                   run->n, wfmt);
 351
 352       tab_double (table, row_headers +i, 5, 0,
 353                   run->runs, &F_8_0);
 354
 355       tab_double (table, row_headers +i, 6, 0,
 356                   z, 0);
 357
 358       tab_double (table, row_headers +i, 7, 0,
 359                   2.0 * gsl_cdf_ugaussian_P (z), 0);
 360     }
 361
 362   switch  ( rt->cp_mode)
 363     {
 364     case CP_CUSTOM:
 365       tab_text (table,  0, column_headers ,
 366                 TAT_TITLE | TAB_LEFT , _("Test Value "));
 367       break;
 368     case CP_MODE:
 369       tab_text (table,  0, column_headers ,
 370                 TAT_TITLE | TAB_LEFT , _("Test Value (mode)"));
 371       break;
 372     case CP_MEAN:
 373       tab_text (table,  0, column_headers ,
 374                 TAT_TITLE | TAB_LEFT , _("Test Value (mean)"));
 375       break;
 376     case CP_MEDIAN:
 377       tab_text (table,  0, column_headers ,
 378                 TAT_TITLE | TAB_LEFT , _("Test Value (median)"));
 379       break;
 380     }
 381
 382   tab_text (table,  0, column_headers + 1,
 383             TAT_TITLE | TAB_LEFT , _("Cases < Test Value"));
 384
 385   tab_text (table,  0, column_headers + 2,
 386             TAT_TITLE | TAB_LEFT , _("Cases >= Test Value"));
 387
 388   tab_text (table,  0, column_headers + 3,
 389             TAT_TITLE | TAB_LEFT , _("Total Cases"));
 390
 391   tab_text (table,  0, column_headers + 4,
 392             TAT_TITLE | TAB_LEFT , _("Number of Runs"));
 393
 394   tab_text (table,  0, column_headers + 5,
 395             TAT_TITLE | TAB_LEFT , _("Z"));
 396
 397   tab_text (table,  0, column_headers + 6,
 398             TAT_TITLE | TAB_LEFT , _("Asymp. Sig. (2-tailed)"));
 399
 400   tab_submit (table);
 401 }
 402
 403