src/language/stats/npar.q

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