GLM: Disable the interactions option until it actually gives the right results
[pspp] / 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