Merge "master" into "output".
[pspp-builds.git] / src / language / stats / npar.q
1 /* PSPP - a program for statistical analysis. -*-c-*-
2    Copyright (C) 2006, 2008, 2009 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 #include <config.h>
18
19 #include <language/stats/npar.h>
20
21 #include <math.h>
22
23 #include <data/case.h>
24 #include <data/casegrouper.h>
25 #include <data/casereader.h>
26 #include <data/dictionary.h>
27 #include <data/procedure.h>
28 #include <language/command.h>
29 #include <language/lexer/lexer.h>
30 #include <language/lexer/variable-parser.h>
31 #include <language/stats/binomial.h>
32 #include <language/stats/chisquare.h>
33 #include <language/stats/wilcoxon.h>
34 #include <language/stats/sign.h>
35 #include <libpspp/cast.h>
36 #include <libpspp/hash.h>
37 #include <libpspp/pool.h>
38 #include <libpspp/taint.h>
39 #include <math/moments.h>
40
41 #include "npar-summary.h"
42
43 #include "gettext.h"
44 #define _(msgid) gettext (msgid)
45
46 /* (headers) */
47
48 /* (specification)
49    "NPAR TESTS" (npar_):
50    +chisquare=custom;
51    +binomial=custom;
52    +wilcoxon=custom;
53    +mcnemar=custom;
54    +sign=custom;
55    +cochran=varlist;
56    +friedman=varlist;
57    +kendall=varlist;
58    missing=miss:!analysis/listwise,
59    incl:include/!exclude;
60    method=custom;
61    +statistics[st_]=descriptives,quartiles,all.
62 */
63 /* (declarations) */
64 /* (functions) */
65
66
67 static struct cmd_npar_tests cmd;
68
69
70 struct npar_specs
71 {
72   struct pool *pool;
73   struct npar_test **test;
74   size_t n_tests;
75
76   const struct variable ** vv; /* Compendium of all variables
77                                   (those mentioned on ANY subcommand */
78   int n_vars; /* Number of variables in vv */
79
80   enum mv_class filter;    /* Missing values to filter. */
81
82   bool descriptives;       /* Descriptive statistics should be calculated */
83   bool quartiles;          /* Quartiles should be calculated */
84
85   bool exact;  /* Whether exact calculations have been requested */
86   double timer;   /* Maximum time (in minutes) to wait for exact calculations */
87 };
88
89 static void one_sample_insert_variables (const struct npar_test *test,
90                                          struct const_hsh_table *variables);
91
92 static void two_sample_insert_variables (const struct npar_test *test,
93                                          struct const_hsh_table *variables);
94
95
96
97 static void
98 npar_execute(struct casereader *input,
99              const struct npar_specs *specs,
100              const struct dataset *ds)
101 {
102   int t;
103   struct descriptives *summary_descriptives = NULL;
104
105   for ( t = 0 ; t < specs->n_tests; ++t )
106     {
107       const struct npar_test *test = specs->test[t];
108       if ( NULL == test->execute )
109         {
110           msg (SW, _("NPAR subcommand not currently implemented."));
111           continue;
112         }
113       test->execute (ds, casereader_clone (input), specs->filter, test, specs->exact, specs->timer);
114     }
115
116   if ( specs->descriptives )
117     {
118       summary_descriptives = xnmalloc (sizeof (*summary_descriptives),
119                                        specs->n_vars);
120
121       npar_summary_calc_descriptives (summary_descriptives,
122                                       casereader_clone (input),
123                                       dataset_dict (ds),
124                                       specs->vv, specs->n_vars,
125                                       specs->filter);
126     }
127
128   if ( (specs->descriptives || specs->quartiles)
129        && !taint_has_tainted_successor (casereader_get_taint (input)) )
130     do_summary_box (summary_descriptives, specs->vv, specs->n_vars );
131
132   free (summary_descriptives);
133   casereader_destroy (input);
134 }
135
136 int
137 cmd_npar_tests (struct lexer *lexer, struct dataset *ds)
138 {
139   bool ok;
140   int i;
141   struct npar_specs npar_specs = {0, 0, 0, 0, 0, 0, 0, 0, 0, 0};
142   struct const_hsh_table *var_hash;
143   struct casegrouper *grouper;
144   struct casereader *input, *group;
145
146   npar_specs.pool = pool_create ();
147
148   var_hash = const_hsh_create_pool (npar_specs.pool, 0,
149                                     compare_vars_by_name, hash_var_by_name,
150                                     NULL, NULL);
151
152   if ( ! parse_npar_tests (lexer, ds, &cmd, &npar_specs) )
153     {
154       pool_destroy (npar_specs.pool);
155       return CMD_FAILURE;
156     }
157
158   for (i = 0; i < npar_specs.n_tests; ++i )
159     {
160       const struct npar_test *test = npar_specs.test[i];
161       test->insert_variables (test, var_hash);
162     }
163
164   npar_specs.vv = (const struct variable **) const_hsh_sort (var_hash);
165   npar_specs.n_vars = const_hsh_count (var_hash);
166
167   if ( cmd.sbc_statistics )
168     {
169       int i;
170
171       for ( i = 0 ; i < NPAR_ST_count; ++i )
172         {
173           if ( cmd.a_statistics[i] )
174             {
175               switch ( i )
176                 {
177                 case NPAR_ST_DESCRIPTIVES:
178                   npar_specs.descriptives = true;
179                   break;
180                 case NPAR_ST_QUARTILES:
181                   npar_specs.quartiles = true;
182                   break;
183                 case NPAR_ST_ALL:
184                   npar_specs.quartiles = true;
185                   npar_specs.descriptives = true;
186                   break;
187                 default:
188                   NOT_REACHED();
189                 };
190             }
191         }
192     }
193
194   npar_specs.filter = cmd.incl == NPAR_EXCLUDE ? MV_ANY : MV_SYSTEM;
195
196   input = proc_open (ds);
197   if ( cmd.miss == NPAR_LISTWISE )
198     {
199       input = casereader_create_filter_missing (input,
200                                                 npar_specs.vv,
201                                                 npar_specs.n_vars,
202                                                 npar_specs.filter,
203                                                 NULL, NULL);
204     }
205
206
207   grouper = casegrouper_create_splits (input, dataset_dict (ds));
208   while (casegrouper_get_next_group (grouper, &group))
209     npar_execute (group, &npar_specs, ds);
210   ok = casegrouper_destroy (grouper);
211   ok = proc_commit (ds) && ok;
212
213   const_hsh_destroy (var_hash);
214
215   pool_destroy (npar_specs.pool);
216
217   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
218 }
219
220 int
221 npar_custom_chisquare (struct lexer *lexer, struct dataset *ds,
222                        struct cmd_npar_tests *cmd UNUSED, void *aux )
223 {
224   struct npar_specs *specs = aux;
225
226   struct chisquare_test *cstp = pool_alloc(specs->pool, sizeof(*cstp));
227   struct one_sample_test *tp = &cstp->parent;
228   struct npar_test *nt = &tp->parent;
229
230   nt->execute = chisquare_execute;
231   nt->insert_variables = one_sample_insert_variables;
232
233   if (!parse_variables_const_pool (lexer, specs->pool, dataset_dict (ds),
234                                    &tp->vars, &tp->n_vars,
235                                    PV_NO_SCRATCH | PV_NO_DUPLICATE))
236     {
237       return 2;
238     }
239
240   cstp->ranged = false;
241
242   if ( lex_match (lexer, '('))
243     {
244       cstp->ranged = true;
245       if ( ! lex_force_num (lexer)) return 0;
246       cstp->lo = lex_integer (lexer);
247       lex_get (lexer);
248       lex_force_match (lexer, ',');
249       if (! lex_force_num (lexer) ) return 0;
250       cstp->hi = lex_integer (lexer);
251       if ( cstp->lo >= cstp->hi )
252         {
253           msg(ME,
254               _("The specified value of HI (%d) is "
255                 "lower than the specified value of LO (%d)"),
256               cstp->hi, cstp->lo);
257           return 0;
258         }
259       lex_get (lexer);
260       if (! lex_force_match (lexer, ')')) return 0;
261     }
262
263   cstp->n_expected = 0;
264   cstp->expected = NULL;
265   if ( lex_match (lexer, '/') )
266     {
267       if ( lex_match_id (lexer, "EXPECTED") )
268         {
269           lex_force_match (lexer, '=');
270           if ( ! lex_match_id (lexer, "EQUAL") )
271             {
272               double f;
273               int n;
274               while ( lex_is_number(lexer) )
275                 {
276                   int i;
277                   n = 1;
278                   f = lex_number (lexer);
279                   lex_get (lexer);
280                   if ( lex_match (lexer, '*'))
281                     {
282                       n = f;
283                       f = lex_number (lexer);
284                       lex_get (lexer);
285                     }
286                   lex_match (lexer, ',');
287
288                   cstp->n_expected += n;
289                   cstp->expected = pool_realloc (specs->pool,
290                                                  cstp->expected,
291                                                  sizeof(double) *
292                                                  cstp->n_expected);
293                   for ( i = cstp->n_expected - n ;
294                         i < cstp->n_expected;
295                         ++i )
296                     cstp->expected[i] = f;
297
298                 }
299             }
300         }
301       else
302         lex_put_back (lexer, '/');
303     }
304
305   if ( cstp->ranged && cstp->n_expected > 0 &&
306        cstp->n_expected != cstp->hi - cstp->lo + 1 )
307     {
308       msg(ME,
309           _("%d expected values were given, but the specified "
310             "range (%d-%d) requires exactly %d values."),
311           cstp->n_expected, cstp->lo, cstp->hi,
312           cstp->hi - cstp->lo +1);
313       return 0;
314     }
315
316   specs->n_tests++;
317   specs->test = pool_realloc (specs->pool,
318                               specs->test,
319                               sizeof(*specs->test) * specs->n_tests);
320
321   specs->test[specs->n_tests - 1] = nt;
322
323   return 1;
324 }
325
326
327 int
328 npar_custom_binomial (struct lexer *lexer, struct dataset *ds,
329                       struct cmd_npar_tests *cmd UNUSED, void *aux)
330 {
331   struct npar_specs *specs = aux;
332   struct binomial_test *btp = pool_alloc(specs->pool, sizeof(*btp));
333   struct one_sample_test *tp = &btp->parent;
334   struct npar_test *nt = &tp->parent;
335
336   nt->execute = binomial_execute;
337   nt->insert_variables = one_sample_insert_variables;
338
339   btp->category1 = btp->category2 = btp->cutpoint = SYSMIS;
340
341   if ( lex_match(lexer, '(') )
342     {
343       if ( lex_force_num (lexer) )
344         {
345           btp->p = lex_number (lexer);
346           lex_get (lexer);
347           lex_force_match (lexer, ')');
348         }
349       else
350         return 0;
351     }
352
353   if ( lex_match (lexer, '=') )
354     {
355       if (parse_variables_const_pool (lexer, specs->pool, dataset_dict (ds),
356                                       &tp->vars, &tp->n_vars,
357                                       PV_NUMERIC | PV_NO_SCRATCH | PV_NO_DUPLICATE) )
358         {
359           if ( lex_match (lexer, '('))
360             {
361               lex_force_num (lexer);
362               btp->category1 = lex_number (lexer);
363               lex_get (lexer);
364               if ( lex_match (lexer, ','))
365                 {
366                   if ( ! lex_force_num (lexer) ) return 2;
367                   btp->category2 = lex_number (lexer);
368                   lex_get (lexer);
369                 }
370               else
371                 {
372                   btp->cutpoint = btp->category1;
373                 }
374
375               lex_force_match (lexer, ')');
376             }
377         }
378       else
379         return 2;
380     }
381
382   specs->n_tests++;
383   specs->test = pool_realloc (specs->pool,
384                               specs->test,
385                               sizeof(*specs->test) * specs->n_tests);
386
387   specs->test[specs->n_tests - 1] = nt;
388
389   return 1;
390 }
391
392
393 bool parse_two_sample_related_test (struct lexer *lexer,
394                                     const struct dictionary *dict,
395                                     struct cmd_npar_tests *cmd,
396                                     struct two_sample_test *test_parameters,
397                                     struct pool *pool
398                                     );
399
400
401 bool
402 parse_two_sample_related_test (struct lexer *lexer,
403                                const struct dictionary *dict,
404                                struct cmd_npar_tests *cmd UNUSED,
405                                struct two_sample_test *test_parameters,
406                                struct pool *pool
407                                )
408 {
409   int n = 0;
410   bool paired = false;
411   bool with = false;
412   const struct variable **vlist1;
413   size_t n_vlist1;
414
415   const struct variable **vlist2;
416   size_t n_vlist2;
417
418   test_parameters->parent.insert_variables = two_sample_insert_variables;
419
420   if (!parse_variables_const_pool (lexer, pool,
421                                    dict,
422                                    &vlist1, &n_vlist1,
423                                    PV_NUMERIC | PV_NO_SCRATCH | PV_NO_DUPLICATE) )
424     return false;
425
426   if ( lex_match(lexer, T_WITH))
427     {
428       with = true;
429       if ( !parse_variables_const_pool (lexer, pool, dict,
430                                         &vlist2, &n_vlist2,
431                                         PV_NUMERIC | PV_NO_SCRATCH | PV_NO_DUPLICATE) )
432         return false;
433
434       paired = (lex_match (lexer, '(') &&
435                 lex_match_id (lexer, "PAIRED") && lex_match (lexer, ')'));
436     }
437
438
439   if ( with )
440     {
441       if (paired)
442         {
443           if ( n_vlist1 != n_vlist2)
444             msg (SE, _("PAIRED was specified but the number of variables "
445                        "preceding WITH (%zu) did not match the number "
446                        "following (%zu)."), n_vlist1, n_vlist2);
447
448           test_parameters->n_pairs = n_vlist1 ;
449         }
450       else
451         {
452           test_parameters->n_pairs = n_vlist1 * n_vlist2;
453         }
454     }
455   else
456     {
457       test_parameters->n_pairs = (n_vlist1 * (n_vlist1 - 1)) / 2 ;
458     }
459
460   test_parameters->pairs =
461     pool_alloc (pool, sizeof ( variable_pair) * test_parameters->n_pairs);
462
463   if ( with )
464     {
465       if (paired)
466         {
467           int i;
468           assert (n_vlist1 == n_vlist2);
469           for ( i = 0 ; i < n_vlist1; ++i )
470             {
471               test_parameters->pairs[n][1] = vlist1[i];
472               test_parameters->pairs[n][0] = vlist2[i];
473               n++;
474             }
475         }
476       else
477         {
478           int i,j;
479           for ( i = 0 ; i < n_vlist1; ++i )
480             {
481               for ( j = 0 ; j < n_vlist2; ++j )
482                 {
483                   test_parameters->pairs[n][1] = vlist1[i];
484                   test_parameters->pairs[n][0] = vlist2[j];
485                   n++;
486                 }
487             }
488         }
489     }
490   else
491     {
492       int i,j;
493       for ( i = 0 ; i < n_vlist1 - 1; ++i )
494         {
495           for ( j = i + 1 ; j < n_vlist1; ++j )
496             {
497               assert ( n < test_parameters->n_pairs);
498               test_parameters->pairs[n][1] = vlist1[i];
499               test_parameters->pairs[n][0] = vlist1[j];
500               n++;
501             }
502         }
503     }
504
505   assert ( n == test_parameters->n_pairs);
506
507   return true;
508 }
509
510 int
511 npar_custom_wilcoxon (struct lexer *lexer,
512                       struct dataset *ds,
513                       struct cmd_npar_tests *cmd, void *aux )
514 {
515   struct npar_specs *specs = aux;
516
517   struct two_sample_test *tp = pool_alloc (specs->pool, sizeof(*tp));
518   struct npar_test *nt = &tp->parent;
519   nt->execute = wilcoxon_execute;
520
521   if (!parse_two_sample_related_test (lexer, dataset_dict (ds), cmd,
522                                       tp, specs->pool) )
523     return 0;
524
525   specs->n_tests++;
526   specs->test = pool_realloc (specs->pool,
527                               specs->test,
528                               sizeof(*specs->test) * specs->n_tests);
529   specs->test[specs->n_tests - 1] = nt;
530
531   return 1;
532 }
533
534 int
535 npar_custom_mcnemar (struct lexer *lexer,
536                      struct dataset *ds,
537                      struct cmd_npar_tests *cmd, void *aux )
538 {
539   struct npar_specs *specs = aux;
540
541   struct two_sample_test *tp = pool_alloc(specs->pool, sizeof(*tp));
542   struct npar_test *nt = &tp->parent;
543   nt->execute = NULL;
544
545
546   if (!parse_two_sample_related_test (lexer, dataset_dict (ds),
547                                       cmd, tp, specs->pool) )
548     return 0;
549
550   specs->n_tests++;
551   specs->test = pool_realloc (specs->pool,
552                               specs->test,
553                               sizeof(*specs->test) * specs->n_tests);
554   specs->test[specs->n_tests - 1] = nt;
555
556   return 1;
557 }
558
559 int
560 npar_custom_sign (struct lexer *lexer, struct dataset *ds,
561                   struct cmd_npar_tests *cmd, void *aux )
562 {
563   struct npar_specs *specs = aux;
564
565   struct two_sample_test *tp = pool_alloc(specs->pool, sizeof(*tp));
566   struct npar_test *nt = &tp->parent;
567
568   nt->execute = sign_execute;
569
570   if (!parse_two_sample_related_test (lexer, dataset_dict (ds), cmd,
571                                       tp, specs->pool) )
572     return 0;
573
574   specs->n_tests++;
575   specs->test = pool_realloc (specs->pool,
576                               specs->test,
577                               sizeof(*specs->test) * specs->n_tests);
578   specs->test[specs->n_tests - 1] = nt;
579
580   return 1;
581 }
582
583 /* Insert the variables for TEST into VAR_HASH */
584 static void
585 one_sample_insert_variables (const struct npar_test *test,
586                              struct const_hsh_table *var_hash)
587 {
588   int i;
589   struct one_sample_test *ost = UP_CAST (test, struct one_sample_test, parent);
590
591   for ( i = 0 ; i < ost->n_vars ; ++i )
592     const_hsh_insert (var_hash, ost->vars[i]);
593 }
594
595 static void
596 two_sample_insert_variables (const struct npar_test *test,
597                              struct const_hsh_table *var_hash)
598 {
599   int i;
600
601   const struct two_sample_test *tst = (const struct two_sample_test *) test;
602
603   for ( i = 0 ; i < tst->n_pairs ; ++i )
604     {
605       variable_pair *pair = &tst->pairs[i];
606
607       const_hsh_insert (var_hash, (*pair)[0]);
608       const_hsh_insert (var_hash, (*pair)[1]);
609     }
610
611 }
612
613
614 static int
615 npar_custom_method (struct lexer *lexer, struct dataset *ds UNUSED,
616                     struct cmd_npar_tests *test UNUSED, void *aux)
617 {
618   struct npar_specs *specs = aux;
619
620   if ( lex_match_id (lexer, "EXACT") )
621     {
622       specs->exact = true;
623       specs->timer = 0.0;
624       if (lex_match_id (lexer, "TIMER"))
625         {
626           specs->timer = 5.0;
627
628           if ( lex_match (lexer, '('))
629             {
630               if ( lex_force_num (lexer) )
631                 {
632                   specs->timer = lex_number (lexer);
633                   lex_get (lexer);
634                 }
635               lex_force_match (lexer, ')');
636             }
637         }
638     }
639
640   return 1;
641 }