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

   1 /* PSPP - a program for statistical analysis. -*-c-*-
   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
  18 #include <config.h>
  19
  20 #include "language/stats/runs.h"
  21
  22 #include <float.h>
  23 #include <gsl/gsl_cdf.h>
  24 #include <math.h>
  25
  26 #include "data/casegrouper.h"
  27 #include "data/casereader.h"
  28 #include "data/casewriter.h"
  29 #include "data/dataset.h"
  30 #include "data/dictionary.h"
  31 #include "data/format.h"
  32 #include "data/subcase.h"
  33 #include "data/variable.h"
  34 #include "libpspp/message.h"
  35 #include "libpspp/misc.h"
  36 #include "math/percentiles.h"
  37 #include "math/sort.h"
  38 #include "output/tab.h"
  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 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 variable `%s'.  "
 169                          "Using %.*g as the threshold value."),
 170                    var_get_name (var), DBL_DIG + 1, run->cutpoint);
 171           }
 172       }
 173       break;
 174     case CP_MEDIAN:
 175       {
 176         for (v = 0; v < otp->n_vars; ++v)
 177           {
 178             double cc = 0.0;
 179             struct ccase *c;
 180             struct run_state *run = &rs[v];
 181             struct casereader *reader = casereader_clone (input);
 182             const struct variable *var = otp->vars[v];
 183             struct casewriter *writer;
 184             struct percentile *median;
 185             struct order_stats *os;
 186             struct subcase sc;
 187             subcase_init_var (&sc, var, SC_ASCEND);
 188             writer = sort_create_writer (&sc, casereader_get_proto (reader));
 189
 190             for (; (c = casereader_read (reader)); )
 191               {
 192                 const union value *val = case_data (c, var);
 193                 const double w = weight ? case_data (c, weight)->f: 1.0;
 194                 if ( var_is_value_missing (var, val, exclude))
 195                   {
 196                     case_unref (c);
 197                     continue;
 198                   }
 199
 200                 cc += w;
 201                 casewriter_write (writer, c);
 202               }
 203             subcase_destroy (&sc);
 204             casereader_destroy (reader);
 205             reader = casewriter_make_reader (writer);
 206
 207             median = percentile_create (0.5, cc);
 208             os = &median->parent;
 209
 210             order_stats_accumulate (&os, 1,
 211                                     reader,
 212                                     weight,
 213                                     var,
 214                                     exclude);
 215
 216             run->cutpoint = percentile_calculate (median, PC_HAVERAGE);
 217             statistic_destroy (&median->parent.parent);
 218           }
 219       }
 220       break;
 221     case CP_MEAN:
 222       {
 223         struct casereader *reader = casereader_clone (input);
 224         for (; (c = casereader_read (reader)); case_unref (c))
 225           {
 226             const double w = weight ? case_data (c, weight)->f: 1.0;
 227             for (v = 0; v < otp->n_vars; ++v)
 228               {
 229                 const struct variable *var = otp->vars[v];
 230                 const union value *val = case_data (c, var);
 231                 const double x = val->f;
 232                 struct run_state *run = &rs[v];
 233
 234                 if ( var_is_value_missing (var, val, exclude))
 235                   continue;
 236
 237                 run->cutpoint += x * w;
 238                 run->n += w;
 239               }
 240           }
 241         casereader_destroy (reader);
 242         for (v = 0; v < otp->n_vars; ++v)
 243           {
 244             struct run_state *run = &rs[v];
 245             run->cutpoint /= run->n;
 246           }
 247       }
 248       break;
 249     case CP_CUSTOM:
 250       {
 251       for (v = 0; v < otp->n_vars; ++v)
 252         {
 253           struct run_state *run = &rs[v];
 254           run->cutpoint = rt->cutpoint;
 255         }
 256       }
 257       break;
 258     }
 259
 260   for (; (c = casereader_read (input)); case_unref (c))
 261     {
 262       const double w = weight ? case_data (c, weight)->f: 1.0;
 263
 264       for (v = 0; v < otp->n_vars; ++v)
 265         {
 266           struct run_state *run = &rs[v];
 267           const struct variable *var = otp->vars[v];
 268           const union value *val = case_data (c, var);
 269           double x = val->f;
 270           double d = x - run->cutpoint;
 271           short sign = 0;
 272
 273           if ( var_is_value_missing (var, val, exclude))
 274             continue;
 275
 276           if (d >= 0)
 277             {
 278               sign = +1;
 279               run->np += w;
 280             }
 281           else
 282             {
 283               sign = -1;
 284               run->nn += w;
 285             }
 286
 287           if (sign != run->last_sign)
 288             run->runs++;
 289
 290           run->last_sign = sign;
 291         }
 292     }
 293   casereader_destroy (input);
 294
 295   for (v = 0; v < otp->n_vars; ++v)
 296     {
 297       struct run_state *run = &rs[v];
 298       run->n = run->np + run->nn;
 299     }
 300
 301   show_runs_result (rt, rs, dict);
 302
 303   free (rs);
 304 }
 305
 306 \f
 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   tab_set_format (table, RC_WEIGHT, wfmt);
 322
 323   tab_headers (table, row_headers, 0, column_headers, 0);
 324
 325   tab_title (table, _("Runs Test"));
 326
 327   /* Box around the table and vertical lines inside*/
 328   tab_box (table, TAL_2, TAL_2, -1, TAL_1,
 329            0,  0, tab_nc (table) - 1, tab_nr (table) - 1 );
 330
 331   tab_hline (table, TAL_2, 0, tab_nc (table) -1, column_headers);
 332   tab_vline (table, TAL_2, row_headers, 0, tab_nr (table) - 1);
 333
 334   for (i = 0 ; i < otp->n_vars; ++i)
 335     {
 336       const struct run_state *run = &rs[i];
 337
 338       double z = runs_statistic (run);
 339
 340       tab_text (table,  row_headers + i, 0,
 341                 TAT_TITLE | TAB_CENTER ,
 342                 var_to_string (otp->vars[i]));
 343
 344       tab_double (table, row_headers +i, 1, 0,
 345                   run->cutpoint, NULL, RC_OTHER);
 346
 347       tab_double (table, row_headers +i, 2, 0,
 348                   run->nn, NULL, RC_WEIGHT);
 349
 350       tab_double (table, row_headers +i, 3, 0,
 351                   run->np, NULL, RC_WEIGHT);
 352
 353       tab_double (table, row_headers +i, 4, 0,
 354                   run->n, NULL, RC_WEIGHT);
 355
 356       tab_double (table, row_headers +i, 5, 0,
 357                   run->runs, NULL, RC_INTEGER);
 358
 359       tab_double (table, row_headers +i, 6, 0,
 360                   z, NULL, RC_OTHER);
 361
 362       tab_double (table, row_headers +i, 7, 0,
 363                   2.0 * (1.0 - gsl_cdf_ugaussian_P (z)), NULL, RC_PVALUE);
 364     }
 365
 366   switch  ( rt->cp_mode)
 367     {
 368     case CP_CUSTOM:
 369       tab_text (table,  0, column_headers ,
 370                 TAT_TITLE | TAB_LEFT , _("Test Value"));
 371       break;
 372     case CP_MODE:
 373       tab_text (table,  0, column_headers ,
 374                 TAT_TITLE | TAB_LEFT , _("Test Value (mode)"));
 375       break;
 376     case CP_MEAN:
 377       tab_text (table,  0, column_headers ,
 378                 TAT_TITLE | TAB_LEFT , _("Test Value (mean)"));
 379       break;
 380     case CP_MEDIAN:
 381       tab_text (table,  0, column_headers ,
 382                 TAT_TITLE | TAB_LEFT , _("Test Value (median)"));
 383       break;
 384     }
 385
 386   tab_text (table,  0, column_headers + 1,
 387             TAT_TITLE | TAB_LEFT , _("Cases < Test Value"));
 388
 389   tab_text (table,  0, column_headers + 2,
 390             TAT_TITLE | TAB_LEFT , _("Cases ≥ Test Value"));
 391
 392   tab_text (table,  0, column_headers + 3,
 393             TAT_TITLE | TAB_LEFT , _("Total Cases"));
 394
 395   tab_text (table,  0, column_headers + 4,
 396             TAT_TITLE | TAB_LEFT , _("Number of Runs"));
 397
 398   tab_text (table,  0, column_headers + 5,
 399             TAT_TITLE | TAB_LEFT , _("Z"));
 400
 401   tab_text (table,  0, column_headers + 6,
 402             TAT_TITLE | TAB_LEFT , _("Asymp. Sig. (2-tailed)"));
 403
 404   tab_submit (table);
 405 }
 406
 407