Cope with new location of upstream files
[pspp] / 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