src/language/stats/t-test.q

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 1997-9, 2000, 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 <gsl/gsl_cdf.h>
  20 #include <math.h>
  21 #include <stdint.h>
  22 #include <stdio.h>
  23 #include <stdlib.h>
  24
  25 #include <data/case.h>
  26 #include <data/casegrouper.h>
  27 #include <data/casereader.h>
  28 #include <data/dictionary.h>
  29 #include <data/procedure.h>
  30 #include <data/value-labels.h>
  31 #include <data/variable.h>
  32 #include <language/command.h>
  33 #include <language/dictionary/split-file.h>
  34 #include <language/lexer/lexer.h>
  35 #include <language/lexer/value-parser.h>
  36 #include <libpspp/array.h>
  37 #include <libpspp/assertion.h>
  38 #include <libpspp/compiler.h>
  39 #include <libpspp/hash.h>
  40 #include <libpspp/message.h>
  41 #include <libpspp/misc.h>
  42 #include <libpspp/str.h>
  43 #include <libpspp/taint.h>
  44 #include <math/group-proc.h>
  45 #include <math/levene.h>
  46 #include <output/manager.h>
  47 #include <output/table.h>
  48 #include <data/format.h>
  49
  50 #include "xalloc.h"
  51 #include "xmemdup0.h"
  52
  53 #include "gettext.h"
  54 #define _(msgid) gettext (msgid)
  55
  56 /* (headers) */
  57
  58 /* (specification)
  59    "T-TEST" (tts_):
  60    +groups=custom;
  61    testval=double;
  62    +variables=varlist("PV_NO_SCRATCH | PV_NUMERIC");
  63    +pairs=custom;
  64    missing=miss:!analysis/listwise,
  65    incl:include/!exclude;
  66    +format=fmt:!labels/nolabels;
  67    criteria=:cin(d:criteria,"%s > 0. && %s < 1.").
  68 */
  69 /* (declarations) */
  70 /* (functions) */
  71
  72 enum comparison
  73   {
  74     CMP_LE,
  75     CMP_EQ,
  76   };
  77
  78 /* A pair of variables to be compared. */
  79 struct pair
  80   {
  81     const struct variable *v[2]; /* The paired variables. */
  82     double n;             /* The number of valid variable pairs */
  83     double sum[2];        /* The sum of the members */
  84     double ssq[2];        /* sum of squares of the members */
  85     double std_dev[2];    /* Std deviation of the members */
  86     double s_std_dev[2];  /* Sample Std deviation of the members */
  87     double mean[2];       /* The means of the members */
  88     double correlation;   /* Correlation coefficient between the variables. */
  89     double sum_of_diffs;  /* The sum of the differences */
  90     double sum_of_prod;   /* The sum of the products */
  91     double mean_diff;     /* The mean of the differences */
  92     double ssq_diffs;     /* The sum of the squares of the differences */
  93     double std_dev_diff;  /* The std deviation of the differences */
  94   };
  95
  96 /* Which mode was T-TEST invoked */
  97 enum t_test_mode {
  98   T_1_SAMPLE,                   /* One-sample tests. */
  99   T_IND_SAMPLES,                /* Independent-sample tests. */
 100   T_PAIRED                      /* Paired-sample tests. */
 101 };
 102
 103 /* Total state of a T-TEST procedure. */
 104 struct t_test_proc
 105   {
 106     enum t_test_mode mode;      /* Mode that T-TEST was invoked in. */
 107     double criteria;            /* Confidence interval in (0, 1). */
 108     enum mv_class exclude;      /* Classes of missing values to exclude. */
 109     bool listwise_missing;      /* Drop whole case if one missing var? */
 110     struct fmt_spec weight_format; /* Format of weight variable. */
 111
 112     /* Dependent variables. */
 113     const struct variable **vars;
 114     size_t n_vars;
 115
 116     /* For mode == T_1_SAMPLE. */
 117     double testval;
 118
 119     /* For mode == T_PAIRED only. */
 120     struct pair *pairs;
 121     size_t n_pairs;
 122
 123     /* For mode == T_IND_SAMPLES only. */
 124     struct variable *indep_var; /* Independent variable. */
 125     enum comparison criterion;  /* Type of comparison. */
 126     double critical_value;      /* CMP_LE only: Grouping threshold value. */
 127     union value g_value[2];     /* CMP_EQ only: Per-group indep var values. */
 128   };
 129
 130 /* Statistics Summary Box */
 131 struct ssbox
 132   {
 133     struct tab_table *t;
 134     void (*populate) (struct ssbox *, struct t_test_proc *);
 135     void (*finalize) (struct ssbox *);
 136   };
 137
 138 static void ssbox_create (struct ssbox *, struct t_test_proc *);
 139 static void ssbox_populate (struct ssbox *, struct t_test_proc *);
 140 static void ssbox_finalize (struct ssbox *);
 141
 142 /* Paired Samples Correlation box */
 143 static void pscbox (struct t_test_proc *);
 144
 145
 146 /* Test Results Box. */
 147 struct trbox {
 148   struct tab_table *t;
 149   void (*populate) (struct trbox *, struct t_test_proc *);
 150   void (*finalize) (struct trbox *);
 151   };
 152
 153 static void trbox_create (struct trbox *, struct t_test_proc *);
 154 static void trbox_populate (struct trbox *, struct t_test_proc *);
 155 static void trbox_finalize (struct trbox *);
 156
 157 static void calculate (struct t_test_proc *, struct casereader *,
 158                        const struct dataset *);
 159
 160 static int compare_group_binary (const struct group_statistics *a,
 161                                  const struct group_statistics *b,
 162                                  const struct t_test_proc *);
 163 static unsigned hash_group_binary (const struct group_statistics *g,
 164                                    const struct t_test_proc *p);
 165
 166 int
 167 cmd_t_test (struct lexer *lexer, struct dataset *ds)
 168 {
 169   struct cmd_t_test cmd;
 170   struct t_test_proc proc;
 171   struct casegrouper *grouper;
 172   struct casereader *group;
 173   struct variable *wv;
 174   bool ok = false;
 175
 176   proc.pairs = NULL;
 177   proc.n_pairs = 0;
 178   proc.vars = NULL;
 179   proc.indep_var = NULL;
 180   if (!parse_t_test (lexer, ds, &cmd, &proc))
 181     goto parse_failed;
 182
 183   wv = dict_get_weight (dataset_dict (ds));
 184   proc.weight_format = wv ? *var_get_print_format (wv) : F_8_0;
 185
 186   if ((cmd.sbc_testval != 0) + (cmd.sbc_groups != 0) + (cmd.sbc_pairs != 0)
 187       != 1)
 188     {
 189       msg (SE, _("Exactly one of TESTVAL, GROUPS and PAIRS subcommands "
 190                  "must be specified."));
 191       goto done;
 192     }
 193
 194   proc.mode = (cmd.sbc_testval ? T_1_SAMPLE
 195                : cmd.sbc_groups ? T_IND_SAMPLES
 196                : T_PAIRED);
 197   proc.criteria = cmd.sbc_criteria ? cmd.criteria : 0.95;
 198   proc.exclude = cmd.incl != TTS_INCLUDE ? MV_ANY : MV_SYSTEM;
 199   proc.listwise_missing = cmd.miss == TTS_LISTWISE;
 200
 201   if (proc.mode == T_1_SAMPLE)
 202     proc.testval = cmd.n_testval[0];
 203
 204   if (proc.mode == T_PAIRED)
 205     {
 206       size_t i, j;
 207
 208       if (cmd.sbc_variables)
 209         {
 210           msg (SE, _("VARIABLES subcommand may not be used with PAIRS."));
 211           goto done;
 212         }
 213
 214       /* Fill proc.vars with the unique variables from pairs. */
 215       proc.n_vars = proc.n_pairs * 2;
 216       proc.vars = xmalloc (sizeof *proc.vars * proc.n_vars);
 217       for (i = j = 0; i < proc.n_pairs; i++)
 218         {
 219           proc.vars[j++] = proc.pairs[i].v[0];
 220           proc.vars[j++] = proc.pairs[i].v[1];
 221         }
 222       proc.n_vars = sort_unique (proc.vars, proc.n_vars, sizeof *proc.vars,
 223                                  compare_var_ptrs_by_name, NULL);
 224     }
 225   else
 226     {
 227       if (!cmd.n_variables)
 228         {
 229           msg (SE, _("One or more VARIABLES must be specified."));
 230           goto done;
 231         }
 232       proc.n_vars = cmd.n_variables;
 233       proc.vars = cmd.v_variables;
 234       cmd.v_variables = NULL;
 235     }
 236
 237   /* Data pass. */
 238   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
 239   while (casegrouper_get_next_group (grouper, &group))
 240     calculate (&proc, group, ds);
 241   ok = casegrouper_destroy (grouper);
 242   ok = proc_commit (ds) && ok;
 243
 244   if (proc.mode == T_IND_SAMPLES)
 245     {
 246       int v;
 247       /* Destroy any group statistics we created */
 248       for (v = 0; v < proc.n_vars; v++)
 249         {
 250           struct group_proc *grpp = group_proc_get (proc.vars[v]);
 251           hsh_destroy (grpp->group_hash);
 252         }
 253     }
 254
 255 done:
 256   free_t_test (&cmd);
 257 parse_failed:
 258   if (proc.indep_var != NULL)
 259     {
 260       int width = var_get_width (proc.indep_var);
 261       value_destroy (&proc.g_value[0], width);
 262       value_destroy (&proc.g_value[1], width);
 263     }
 264   free (proc.vars);
 265   free (proc.pairs);
 266   return ok ? CMD_SUCCESS : CMD_FAILURE;
 267 }
 268
 269 static int
 270 tts_custom_groups (struct lexer *lexer, struct dataset *ds,
 271                    struct cmd_t_test *cmd UNUSED, void *proc_)
 272 {
 273   struct t_test_proc *proc = proc_;
 274   int n_values;
 275   int width;
 276
 277   lex_match (lexer, '=');
 278
 279   proc->indep_var = parse_variable (lexer, dataset_dict (ds));
 280   if (proc->indep_var == NULL)
 281     {
 282       lex_error (lexer, "expecting variable name in GROUPS subcommand");
 283       return 0;
 284     }
 285   width = var_get_width (proc->indep_var);
 286   value_init (&proc->g_value[0], width);
 287   value_init (&proc->g_value[1], width);
 288
 289   if (!lex_match (lexer, '('))
 290     n_values = 0;
 291   else
 292     {
 293       if (!parse_value (lexer, &proc->g_value[0], width))
 294         return 0;
 295       lex_match (lexer, ',');
 296       if (lex_match (lexer, ')'))
 297         n_values = 1;
 298       else
 299         {
 300           if (!parse_value (lexer, &proc->g_value[1], width)
 301               || !lex_force_match (lexer, ')'))
 302             return 0;
 303           n_values = 2;
 304         }
 305     }
 306
 307   if (var_is_numeric (proc->indep_var))
 308     {
 309       proc->criterion = n_values == 1 ? CMP_LE : CMP_EQ;
 310       if (n_values == 1)
 311         proc->critical_value = proc->g_value[0].f;
 312       else if (n_values == 0)
 313         {
 314           proc->g_value[0].f = 1;
 315           proc->g_value[1].f = 2;
 316         }
 317     }
 318   else
 319     {
 320       proc->criterion = CMP_EQ;
 321       if (n_values != 2)
 322         {
 323           msg (SE, _("When applying GROUPS to a string variable, two "
 324                      "values must be specified."));
 325           return 0;
 326         }
 327     }
 328   return 1;
 329 }
 330
 331 static void
 332 add_pair (struct t_test_proc *proc,
 333           const struct variable *v0, const struct variable *v1)
 334 {
 335   struct pair *p = &proc->pairs[proc->n_pairs++];
 336   p->v[0] = v0;
 337   p->v[1] = v1;
 338 }
 339
 340 static int
 341 tts_custom_pairs (struct lexer *lexer, struct dataset *ds,
 342                   struct cmd_t_test *cmd UNUSED, void *proc_)
 343 {
 344   struct t_test_proc *proc = proc_;
 345
 346   const struct variable **vars1 = NULL;
 347   size_t n_vars1 = 0;
 348
 349   const struct variable **vars2 = NULL;
 350   size_t n_vars2 = 0;
 351
 352   bool paired = false;
 353
 354   size_t n_total_pairs;
 355   size_t i, j;
 356
 357   lex_match (lexer, '=');
 358
 359   if (!parse_variables_const (lexer, dataset_dict (ds), &vars1, &n_vars1,
 360                               PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH))
 361     return 0;
 362
 363   if (lex_match (lexer, T_WITH))
 364     {
 365       if (!parse_variables_const (lexer, dataset_dict (ds), &vars2, &n_vars2,
 366                                   PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH))
 367         {
 368           free (vars1);
 369           return 0;
 370         }
 371
 372       if (lex_match (lexer, '(')
 373           && lex_match_id (lexer, "PAIRED")
 374           && lex_match (lexer, ')'))
 375         {
 376           paired = true;
 377           if (n_vars1 != n_vars2)
 378             {
 379               msg (SE, _("PAIRED was specified but the number of variables "
 380                          "preceding WITH (%zu) did not match the number "
 381                          "following (%zu)."),
 382                    n_vars1, n_vars2);
 383               free (vars1);
 384               free (vars2);
 385               return 0;
 386             }
 387         }
 388     }
 389   else
 390     {
 391       if (n_vars1 < 2)
 392         {
 393           free (vars1);
 394           msg (SE, _("At least two variables must be specified on PAIRS."));
 395           return 0;
 396         }
 397     }
 398
 399   /* Allocate storage for the new pairs. */
 400   n_total_pairs = proc->n_pairs + (paired ? n_vars1
 401                                    : n_vars2 > 0 ? n_vars1 * n_vars2
 402                                    : n_vars1 * (n_vars1 - 1) / 2);
 403   proc->pairs = xnrealloc (proc->pairs, n_total_pairs, sizeof *proc->pairs);
 404
 405   /* Populate the pairs with the appropriate variables. */
 406   if (paired)
 407     for (i = 0; i < n_vars1; i++)
 408       add_pair (proc, vars1[i], vars2[i]);
 409   else if (n_vars2 > 0)
 410     for (i = 0; i < n_vars1; i++)
 411       for (j = 0; j < n_vars2; j++)
 412         add_pair (proc, vars1[i], vars2[j]);
 413   else
 414     for (i = 0; i < n_vars1; i++)
 415       for (j = i + 1; j < n_vars1; j++)
 416         add_pair (proc, vars1[i], vars1[j]);
 417   assert (proc->n_pairs == n_total_pairs);
 418
 419   free (vars1);
 420   free (vars2);
 421   return 1;
 422 }
 423 \f
 424 /* Implementation of the SSBOX object. */
 425
 426 static void ssbox_base_init (struct ssbox *, int cols, int rows);
 427 static void ssbox_base_finalize (struct ssbox *);
 428 static void ssbox_one_sample_init (struct ssbox *, struct t_test_proc *);
 429 static void ssbox_independent_samples_init (struct ssbox *, struct t_test_proc *);
 430 static void ssbox_paired_init (struct ssbox *, struct t_test_proc *);
 431
 432 /* Factory to create an ssbox. */
 433 static void
 434 ssbox_create (struct ssbox *ssb, struct t_test_proc *proc)
 435 {
 436   switch (proc->mode)
 437     {
 438     case T_1_SAMPLE:
 439       ssbox_one_sample_init (ssb, proc);
 440       break;
 441     case T_IND_SAMPLES:
 442       ssbox_independent_samples_init (ssb, proc);
 443       break;
 444     case T_PAIRED:
 445       ssbox_paired_init (ssb, proc);
 446       break;
 447     default:
 448       NOT_REACHED ();
 449     }
 450 }
 451
 452 /* Despatcher for the populate method */
 453 static void
 454 ssbox_populate (struct ssbox *ssb, struct t_test_proc *proc)
 455 {
 456   ssb->populate (ssb, proc);
 457 }
 458
 459 /* Despatcher for finalize */
 460 static void
 461 ssbox_finalize (struct ssbox *ssb)
 462 {
 463   ssb->finalize (ssb);
 464 }
 465
 466 /* Submit the box and clear up */
 467 static void
 468 ssbox_base_finalize (struct ssbox *ssb)
 469 {
 470   tab_submit (ssb->t);
 471 }
 472
 473 /* Initialize a ssbox struct */
 474 static void
 475 ssbox_base_init (struct ssbox *this, int cols, int rows)
 476 {
 477   this->finalize = ssbox_base_finalize;
 478   this->t = tab_create (cols, rows, 0);
 479
 480   tab_columns (this->t, SOM_COL_DOWN, 1);
 481   tab_headers (this->t, 0, 0, 1, 0);
 482   tab_box (this->t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols - 1, rows - 1);
 483   tab_hline (this->t, TAL_2, 0, cols- 1, 1);
 484   tab_dim (this->t, tab_natural_dimensions, NULL);
 485 }
 486 \f
 487 /* ssbox implementations. */
 488
 489 static void ssbox_one_sample_populate (struct ssbox *, struct t_test_proc *);
 490 static void ssbox_independent_samples_populate (struct ssbox *,
 491                                                 struct t_test_proc *);
 492 static void ssbox_paired_populate (struct ssbox *, struct t_test_proc *);
 493
 494 /* Initialize the one_sample ssbox */
 495 static void
 496 ssbox_one_sample_init (struct ssbox *this, struct t_test_proc *proc)
 497 {
 498   const int hsize = 5;
 499   const int vsize = proc->n_vars + 1;
 500
 501   this->populate = ssbox_one_sample_populate;
 502
 503   ssbox_base_init (this, hsize, vsize);
 504   tab_title (this->t, _("One-Sample Statistics"));
 505   tab_vline (this->t, TAL_2, 1, 0, vsize - 1);
 506   tab_text (this->t, 1, 0, TAB_CENTER | TAT_TITLE, _("N"));
 507   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
 508   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 509   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _("SE. Mean"));
 510 }
 511
 512 /* Initialize the independent samples ssbox */
 513 static void
 514 ssbox_independent_samples_init (struct ssbox *this, struct t_test_proc *proc)
 515 {
 516   int hsize=6;
 517   int vsize = proc->n_vars * 2 + 1;
 518
 519   this->populate = ssbox_independent_samples_populate;
 520
 521   ssbox_base_init (this, hsize, vsize);
 522   tab_vline (this->t, TAL_GAP, 1, 0, vsize - 1);
 523   tab_title (this->t, _("Group Statistics"));
 524   tab_text (this->t, 1, 0, TAB_CENTER | TAT_TITLE,
 525             var_get_name (proc->indep_var));
 526   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _("N"));
 527   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
 528   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 529   tab_text (this->t, 5, 0, TAB_CENTER | TAT_TITLE, _("SE. Mean"));
 530 }
 531
 532 /* Populate the ssbox for independent samples */
 533 static void
 534 ssbox_independent_samples_populate (struct ssbox *ssb,
 535                                     struct t_test_proc *proc)
 536 {
 537   int i;
 538
 539   char *val_lab[2];
 540   double indep_value[2];
 541
 542   char prefix[2][3];
 543
 544   for (i = 0; i < 2; i++)
 545     {
 546       union value *value = &proc->g_value[i];
 547       int width = var_get_width (proc->indep_var);
 548
 549       indep_value[i] = (proc->criterion == CMP_LE ? proc->critical_value
 550                         : value->f);
 551
 552       if (val_type_from_width (width) == VAL_NUMERIC)
 553         {
 554           const char *s = var_lookup_value_label (proc->indep_var, value);
 555           val_lab[i] = s ? xstrdup (s) : xasprintf ("%g", indep_value[i]);
 556         }
 557       else
 558         val_lab[i] = xmemdup0 (value_str (value, width), width);
 559     }
 560
 561   if (proc->criterion == CMP_LE)
 562     {
 563       strcpy (prefix[0], ">=");
 564       strcpy (prefix[1], "<");
 565     }
 566   else
 567     {
 568       strcpy (prefix[0], "");
 569       strcpy (prefix[1], "");
 570     }
 571
 572   for (i = 0; i < proc->n_vars; i++)
 573     {
 574       const struct variable *var = proc->vars[i];
 575       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
 576       int count=0;
 577
 578       tab_text (ssb->t, 0, i * 2 + 1, TAB_LEFT,
 579                 var_get_name (proc->vars[i]));
 580       tab_text_format (ssb->t, 1, i * 2 + 1, TAB_LEFT,
 581                        "%s%s", prefix[0], val_lab[0]);
 582       tab_text_format (ssb->t, 1, i * 2 + 1+ 1, TAB_LEFT,
 583                        "%s%s", prefix[1], val_lab[1]);
 584
 585       /* Fill in the group statistics */
 586       for (count = 0; count < 2; count++)
 587         {
 588           union value search_val;
 589           struct group_statistics *gs;
 590
 591           if (proc->criterion == CMP_LE)
 592             search_val.f = proc->critical_value + (count == 0 ? 1.0 : -1.0);
 593           else
 594             search_val = proc->g_value[count];
 595
 596           gs = hsh_find (grp_hash, &search_val);
 597           assert (gs);
 598
 599           tab_double (ssb->t, 2, i * 2 + count+ 1, TAB_RIGHT, gs->n,
 600                       &proc->weight_format);
 601           tab_double (ssb->t, 3, i * 2 + count+ 1, TAB_RIGHT, gs->mean, NULL);
 602           tab_double (ssb->t, 4, i * 2 + count+ 1, TAB_RIGHT, gs->std_dev,
 603                       NULL);
 604           tab_double (ssb->t, 5, i * 2 + count+ 1, TAB_RIGHT, gs->se_mean,
 605                       NULL);
 606         }
 607     }
 608   free (val_lab[0]);
 609   free (val_lab[1]);
 610 }
 611
 612 /* Initialize the paired values ssbox */
 613 static void
 614 ssbox_paired_init (struct ssbox *this, struct t_test_proc *proc)
 615 {
 616   int hsize = 6;
 617   int vsize = proc->n_pairs * 2 + 1;
 618
 619   this->populate = ssbox_paired_populate;
 620
 621   ssbox_base_init (this, hsize, vsize);
 622   tab_title (this->t, _("Paired Sample Statistics"));
 623   tab_vline (this->t, TAL_GAP, 1, 0, vsize - 1);
 624   tab_vline (this->t, TAL_2, 2, 0, vsize - 1);
 625   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _("Mean"));
 626   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _("N"));
 627   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 628   tab_text (this->t, 5, 0, TAB_CENTER | TAT_TITLE, _("SE. Mean"));
 629 }
 630
 631 /* Populate the ssbox for paired values */
 632 static void
 633 ssbox_paired_populate (struct ssbox *ssb, struct t_test_proc *proc)
 634 {
 635   int i;
 636
 637   for (i = 0; i < proc->n_pairs; i++)
 638     {
 639       struct pair *p = &proc->pairs[i];
 640       int j;
 641
 642       tab_text_format (ssb->t, 0, i * 2 + 1, TAB_LEFT, _("Pair %d"), i);
 643       for (j=0; j < 2; j++)
 644         {
 645           /* Titles */
 646           tab_text (ssb->t, 1, i * 2 + j + 1, TAB_LEFT,
 647                     var_get_name (p->v[j]));
 648
 649           /* Values */
 650           tab_double (ssb->t, 2, i * 2 + j + 1, TAB_RIGHT, p->mean[j], NULL);
 651           tab_double (ssb->t, 3, i * 2 + j + 1, TAB_RIGHT, p->n,
 652                       &proc->weight_format);
 653           tab_double (ssb->t, 4, i * 2 + j + 1, TAB_RIGHT, p->std_dev[j],
 654                       NULL);
 655           tab_double (ssb->t, 5, i * 2 + j + 1, TAB_RIGHT,
 656                      p->std_dev[j] /sqrt (p->n), NULL);
 657         }
 658     }
 659 }
 660
 661 /* Populate the one sample ssbox */
 662 static void
 663 ssbox_one_sample_populate (struct ssbox *ssb, struct t_test_proc *proc)
 664 {
 665   int i;
 666
 667   for (i = 0; i < proc->n_vars; i++)
 668     {
 669       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
 670
 671       tab_text (ssb->t, 0, i + 1, TAB_LEFT, var_get_name (proc->vars[i]));
 672       tab_double (ssb->t, 1, i + 1, TAB_RIGHT, gs->n, &proc->weight_format);
 673       tab_double (ssb->t, 2, i + 1, TAB_RIGHT, gs->mean, NULL);
 674       tab_double (ssb->t, 3, i + 1, TAB_RIGHT, gs->std_dev, NULL);
 675       tab_double (ssb->t, 4, i + 1, TAB_RIGHT, gs->se_mean, NULL);
 676     }
 677 }
 678 \f
 679 /* Implementation of the Test Results box struct */
 680
 681 static void trbox_base_init (struct trbox *, size_t n_vars, int cols);
 682 static void trbox_base_finalize (struct trbox *);
 683 static void trbox_independent_samples_init (struct trbox *,
 684                                             struct t_test_proc *);
 685 static void trbox_independent_samples_populate (struct trbox *,
 686                                                 struct t_test_proc *);
 687 static void trbox_one_sample_init (struct trbox *, struct t_test_proc *);
 688 static void trbox_one_sample_populate (struct trbox *, struct t_test_proc *);
 689 static void trbox_paired_init (struct trbox *, struct t_test_proc *);
 690 static void trbox_paired_populate (struct trbox *, struct t_test_proc *);
 691
 692 /* Create a trbox according to mode*/
 693 static void
 694 trbox_create (struct trbox *trb, struct t_test_proc *proc)
 695 {
 696   switch (proc->mode)
 697     {
 698     case T_1_SAMPLE:
 699       trbox_one_sample_init (trb, proc);
 700       break;
 701     case T_IND_SAMPLES:
 702       trbox_independent_samples_init (trb, proc);
 703       break;
 704     case T_PAIRED:
 705       trbox_paired_init (trb, proc);
 706       break;
 707     default:
 708       NOT_REACHED ();
 709     }
 710 }
 711
 712 /* Populate a trbox according to proc */
 713 static void
 714 trbox_populate (struct trbox *trb, struct t_test_proc *proc)
 715 {
 716   trb->populate (trb, proc);
 717 }
 718
 719 /* Submit and destroy a trbox */
 720 static void
 721 trbox_finalize (struct trbox *trb)
 722 {
 723   trb->finalize (trb);
 724 }
 725
 726 /* Initialize the independent samples trbox */
 727 static void
 728 trbox_independent_samples_init (struct trbox *self,
 729                                 struct t_test_proc *proc)
 730 {
 731   const int hsize = 11;
 732   const int vsize = proc->n_vars * 2 + 3;
 733
 734   assert (self);
 735   self->populate = trbox_independent_samples_populate;
 736
 737   trbox_base_init (self, proc->n_vars * 2, hsize);
 738   tab_title (self->t, _("Independent Samples Test"));
 739   tab_hline (self->t, TAL_1, 2, hsize - 1, 1);
 740   tab_vline (self->t, TAL_2, 2, 0, vsize - 1);
 741   tab_vline (self->t, TAL_1, 4, 0, vsize - 1);
 742   tab_box (self->t, -1, -1, -1, TAL_1, 2, 1, hsize - 2, vsize - 1);
 743   tab_hline (self->t, TAL_1, hsize - 2, hsize - 1, 2);
 744   tab_box (self->t, -1, -1, -1, TAL_1, hsize - 2, 2, hsize - 1, vsize - 1);
 745   tab_joint_text (self->t, 2, 0, 3, 0,
 746                   TAB_CENTER, _("Levene's Test for Equality of Variances"));
 747   tab_joint_text (self->t, 4, 0, hsize- 1, 0,
 748                   TAB_CENTER, _("t-test for Equality of Means"));
 749
 750   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _("F"));
 751   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _("Sig."));
 752   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _("t"));
 753   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _("df"));
 754   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
 755   tab_text (self->t, 7, 2, TAB_CENTER | TAT_TITLE, _("Mean Difference"));
 756   tab_text (self->t, 8, 2, TAB_CENTER | TAT_TITLE, _("Std. Error Difference"));
 757   tab_text (self->t, 9, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
 758   tab_text (self->t, 10, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
 759
 760   tab_joint_text_format (self->t, 9, 1, 10, 1, TAB_CENTER,
 761                          _("%g%% Confidence Interval of the Difference"),
 762                          proc->criteria * 100.0);
 763 }
 764
 765 /* Populate the independent samples trbox */
 766 static void
 767 trbox_independent_samples_populate (struct trbox *self,
 768                                     struct t_test_proc *proc)
 769 {
 770   int i;
 771
 772   for (i = 0; i < proc->n_vars; i++)
 773     {
 774       double p, q;
 775
 776       double t;
 777       double df;
 778
 779       double df1, df2;
 780
 781       double pooled_variance;
 782       double std_err_diff;
 783       double mean_diff;
 784
 785       double se2;
 786
 787       const struct variable *var = proc->vars[i];
 788       struct group_proc *grp_data = group_proc_get (var);
 789
 790       struct hsh_table *grp_hash = grp_data->group_hash;
 791
 792       struct group_statistics *gs0;
 793       struct group_statistics *gs1;
 794
 795       union value search_val;
 796
 797       if (proc->criterion == CMP_LE)
 798         search_val.f = proc->critical_value - 1.0;
 799       else
 800         search_val = proc->g_value[0];
 801
 802       gs0 = hsh_find (grp_hash, &search_val);
 803       assert (gs0);
 804
 805       if (proc->criterion == CMP_LE)
 806         search_val.f = proc->critical_value + 1.0;
 807       else
 808         search_val = proc->g_value[1];
 809
 810       gs1 = hsh_find (grp_hash, &search_val);
 811       assert (gs1);
 812
 813
 814       tab_text (self->t, 0, i * 2 + 3, TAB_LEFT, var_get_name (proc->vars[i]));
 815       tab_text (self->t, 1, i * 2 + 3, TAB_LEFT, _("Equal variances assumed"));
 816       tab_double (self->t, 2, i * 2 + 3, TAB_CENTER, grp_data->levene, NULL);
 817
 818       /* Now work out the significance of the Levene test */
 819       df1 = 1;
 820       df2 = grp_data->ugs.n - 2;
 821       q = gsl_cdf_fdist_Q (grp_data->levene, df1, df2);
 822       tab_double (self->t, 3, i * 2 + 3, TAB_CENTER, q, NULL);
 823
 824       df = gs0->n + gs1->n - 2.0;
 825       tab_double (self->t, 5, i * 2 + 3, TAB_RIGHT, df, NULL);
 826
 827       pooled_variance = (gs0->n * pow2 (gs0->s_std_dev)
 828                          + gs1->n *pow2 (gs1->s_std_dev)) / df ;
 829
 830       t = (gs0->mean - gs1->mean) / sqrt (pooled_variance);
 831       t /= sqrt ((gs0->n + gs1->n) / (gs0->n * gs1->n));
 832
 833       tab_double (self->t, 4, i * 2 + 3, TAB_RIGHT, t, NULL);
 834
 835       p = gsl_cdf_tdist_P (t, df);
 836       q = gsl_cdf_tdist_Q (t, df);
 837
 838       tab_double (self->t, 6, i * 2 + 3, TAB_RIGHT, 2.0 * (t > 0 ? q : p),
 839                   NULL);
 840
 841       mean_diff = gs0->mean - gs1->mean;
 842       tab_double (self->t, 7, i * 2 + 3, TAB_RIGHT, mean_diff, NULL);
 843
 844
 845       std_err_diff = sqrt (pow2 (gs0->se_mean) + pow2 (gs1->se_mean));
 846       tab_double (self->t, 8, i * 2 + 3, TAB_RIGHT, std_err_diff, NULL);
 847
 848       /* Now work out the confidence interval */
 849       q = (1 - proc->criteria)/2.0;  /* 2-tailed test */
 850
 851       t = gsl_cdf_tdist_Qinv (q, df);
 852       tab_double (self->t, 9, i * 2 + 3, TAB_RIGHT,
 853                  mean_diff - t * std_err_diff, NULL);
 854
 855       tab_double (self->t, 10, i * 2 + 3, TAB_RIGHT,
 856                  mean_diff + t * std_err_diff, NULL);
 857
 858
 859       /* Now for the \sigma_1 != \sigma_2 case */
 860       tab_text (self->t, 1, i * 2 + 3 + 1,
 861                 TAB_LEFT, _("Equal variances not assumed"));
 862
 863       se2 = ((pow2 (gs0->s_std_dev) / (gs0->n - 1)) +
 864              (pow2 (gs1->s_std_dev) / (gs1->n - 1)));
 865
 866       t = mean_diff / sqrt (se2);
 867       tab_double (self->t, 4, i * 2 + 3 + 1, TAB_RIGHT, t, NULL);
 868
 869       df = pow2 (se2) / ((pow2 (pow2 (gs0->s_std_dev) / (gs0->n - 1))
 870                           / (gs0->n - 1))
 871                          + (pow2 (pow2 (gs1->s_std_dev) / (gs1->n - 1))
 872                             / (gs1->n - 1)));
 873       tab_double (self->t, 5, i * 2 + 3 + 1, TAB_RIGHT, df, NULL);
 874
 875       p = gsl_cdf_tdist_P (t, df);
 876       q = gsl_cdf_tdist_Q (t, df);
 877
 878       tab_double (self->t, 6, i * 2 + 3 + 1, TAB_RIGHT, 2.0 * (t > 0 ? q : p),
 879                   NULL);
 880
 881       /* Now work out the confidence interval */
 882       q = (1 - proc->criteria) / 2.0;  /* 2-tailed test */
 883
 884       t = gsl_cdf_tdist_Qinv (q, df);
 885
 886       tab_double (self->t, 7, i * 2 + 3 + 1, TAB_RIGHT, mean_diff, NULL);
 887       tab_double (self->t, 8, i * 2 + 3 + 1, TAB_RIGHT, std_err_diff, NULL);
 888       tab_double (self->t, 9, i * 2 + 3 + 1, TAB_RIGHT,
 889                  mean_diff - t * std_err_diff, NULL);
 890       tab_double (self->t, 10, i * 2 + 3 + 1, TAB_RIGHT,
 891                  mean_diff + t * std_err_diff, NULL);
 892     }
 893 }
 894
 895 /* Initialize the paired samples trbox */
 896 static void
 897 trbox_paired_init (struct trbox *self, struct t_test_proc *proc)
 898 {
 899   const int hsize=10;
 900   const int vsize=proc->n_pairs+ 3;
 901
 902   self->populate = trbox_paired_populate;
 903
 904   trbox_base_init (self, proc->n_pairs, hsize);
 905   tab_title (self->t, _("Paired Samples Test"));
 906   tab_hline (self->t, TAL_1, 2, 6, 1);
 907   tab_vline (self->t, TAL_2, 2, 0, vsize - 1);
 908   tab_joint_text (self->t, 2, 0, 6, 0, TAB_CENTER, _("Paired Differences"));
 909   tab_box (self->t, -1, -1, -1, TAL_1, 2, 1, 6, vsize - 1);
 910   tab_box (self->t, -1, -1, -1, TAL_1, 6, 0, hsize - 1, vsize - 1);
 911   tab_hline (self->t, TAL_1, 5, 6, 2);
 912   tab_vline (self->t, TAL_GAP, 6, 0, 1);
 913
 914   tab_joint_text_format (self->t, 5, 1, 6, 1, TAB_CENTER,
 915                          _("%g%% Confidence Interval of the Difference"),
 916                          proc->criteria*100.0);
 917
 918   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _("Mean"));
 919   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
 920   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _("Std. Error Mean"));
 921   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
 922   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
 923   tab_text (self->t, 7, 2, TAB_CENTER | TAT_TITLE, _("t"));
 924   tab_text (self->t, 8, 2, TAB_CENTER | TAT_TITLE, _("df"));
 925   tab_text (self->t, 9, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
 926 }
 927
 928 /* Populate the paired samples trbox */
 929 static void
 930 trbox_paired_populate (struct trbox *trb,
 931                        struct t_test_proc *proc)
 932 {
 933   int i;
 934
 935   for (i = 0; i < proc->n_pairs; i++)
 936     {
 937       struct pair *pair = &proc->pairs[i];
 938       double p, q;
 939       double se_mean;
 940
 941       double n = pair->n;
 942       double t;
 943       double df = n - 1;
 944
 945       tab_text_format (trb->t, 0, i + 3, TAB_LEFT, _("Pair %d"), i);
 946       tab_text_format (trb->t, 1, i + 3, TAB_LEFT, "%s - %s",
 947                        var_get_name (pair->v[0]),
 948                        var_get_name (pair->v[1]));
 949       tab_double (trb->t, 2, i + 3, TAB_RIGHT, pair->mean_diff, NULL);
 950       tab_double (trb->t, 3, i + 3, TAB_RIGHT, pair->std_dev_diff, NULL);
 951
 952       /* SE Mean */
 953       se_mean = pair->std_dev_diff / sqrt (n);
 954       tab_double (trb->t, 4, i + 3, TAB_RIGHT, se_mean, NULL);
 955
 956       /* Now work out the confidence interval */
 957       q = (1 - proc->criteria) / 2.0;  /* 2-tailed test */
 958
 959       t = gsl_cdf_tdist_Qinv (q, df);
 960
 961       tab_double (trb->t, 5, i + 3, TAB_RIGHT,
 962                  pair->mean_diff - t * se_mean, NULL);
 963       tab_double (trb->t, 6, i + 3, TAB_RIGHT,
 964                  pair->mean_diff + t * se_mean, NULL);
 965
 966       t = ((pair->mean[0] - pair->mean[1])
 967            / sqrt ((pow2 (pair->s_std_dev[0]) + pow2 (pair->s_std_dev[1])
 968                     - (2 * pair->correlation
 969                        * pair->s_std_dev[0] * pair->s_std_dev[1]))
 970                    / (n - 1)));
 971
 972       tab_double (trb->t, 7, i + 3, TAB_RIGHT, t, NULL);
 973
 974       /* Degrees of freedom */
 975       tab_double (trb->t, 8, i + 3, TAB_RIGHT, df, &proc->weight_format);
 976
 977       p = gsl_cdf_tdist_P (t, df);
 978       q = gsl_cdf_tdist_P (t, df);
 979
 980       tab_double (trb->t, 9, i + 3, TAB_RIGHT, 2.0 * (t > 0 ? q : p), NULL);
 981     }
 982 }
 983
 984 /* Initialize the one sample trbox */
 985 static void
 986 trbox_one_sample_init (struct trbox *self, struct t_test_proc *proc)
 987 {
 988   const int hsize = 7;
 989   const int vsize = proc->n_vars + 3;
 990
 991   self->populate = trbox_one_sample_populate;
 992
 993   trbox_base_init (self, proc->n_vars, hsize);
 994   tab_title (self->t, _("One-Sample Test"));
 995   tab_hline (self->t, TAL_1, 1, hsize - 1, 1);
 996   tab_vline (self->t, TAL_2, 1, 0, vsize - 1);
 997
 998   tab_joint_text_format (self->t, 1, 0, hsize - 1, 0, TAB_CENTER,
 999                          _("Test Value = %f"), proc->testval);
1000
1001   tab_box (self->t, -1, -1, -1, TAL_1, 1, 1, hsize - 1, vsize - 1);
1002
1003
1004   tab_joint_text_format (self->t, 5, 1, 6, 1, TAB_CENTER,
1005                          _("%g%% Confidence Interval of the Difference"),
1006                          proc->criteria * 100.0);
1007
1008   tab_vline (self->t, TAL_GAP, 6, 1, 1);
1009   tab_hline (self->t, TAL_1, 5, 6, 2);
1010   tab_text (self->t, 1, 2, TAB_CENTER | TAT_TITLE, _("t"));
1011   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _("df"));
1012   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
1013   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _("Mean Difference"));
1014   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _("Lower"));
1015   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _("Upper"));
1016 }
1017
1018 /* Populate the one sample trbox */
1019 static void
1020 trbox_one_sample_populate (struct trbox *trb, struct t_test_proc *proc)
1021 {
1022   int i;
1023
1024   assert (trb->t);
1025
1026   for (i = 0; i < proc->n_vars; i++)
1027     {
1028       double t;
1029       double p, q;
1030       double df;
1031       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
1032
1033       tab_text (trb->t, 0, i + 3, TAB_LEFT, var_get_name (proc->vars[i]));
1034
1035       t = (gs->mean - proc->testval) * sqrt (gs->n) / gs->std_dev;
1036
1037       tab_double (trb->t, 1, i + 3, TAB_RIGHT, t, NULL);
1038
1039       /* degrees of freedom */
1040       df = gs->n - 1;
1041
1042       tab_double (trb->t, 2, i + 3, TAB_RIGHT, df, &proc->weight_format);
1043
1044       p = gsl_cdf_tdist_P (t, df);
1045       q = gsl_cdf_tdist_Q (t, df);
1046
1047       /* Multiply by 2 to get 2-tailed significance, makeing sure we've got
1048          the correct tail*/
1049       tab_double (trb->t, 3, i + 3, TAB_RIGHT, 2.0 * (t > 0 ? q : p), NULL);
1050       tab_double (trb->t, 4, i + 3, TAB_RIGHT, gs->mean_diff, NULL);
1051
1052
1053       q = (1 - proc->criteria) / 2.0;  /* 2-tailed test */
1054       t = gsl_cdf_tdist_Qinv (q, df);
1055
1056       tab_double (trb->t, 5, i + 3, TAB_RIGHT,
1057                  gs->mean_diff - t * gs->se_mean, NULL);
1058       tab_double (trb->t, 6, i + 3, TAB_RIGHT,
1059                  gs->mean_diff + t * gs->se_mean, NULL);
1060     }
1061 }
1062
1063 /* Base initializer for the generalized trbox */
1064 static void
1065 trbox_base_init (struct trbox *self, size_t data_rows, int cols)
1066 {
1067   const size_t rows = 3 + data_rows;
1068
1069   self->finalize = trbox_base_finalize;
1070   self->t = tab_create (cols, rows, 0);
1071   tab_headers (self->t, 0, 0, 3, 0);
1072   tab_box (self->t, TAL_2, TAL_2, TAL_0, TAL_0, 0, 0, cols - 1, rows - 1);
1073   tab_hline (self->t, TAL_2, 0, cols- 1, 3);
1074   tab_dim (self->t, tab_natural_dimensions, NULL);
1075 }
1076
1077 /* Base finalizer for the trbox */
1078 static void
1079 trbox_base_finalize (struct trbox *trb)
1080 {
1081   tab_submit (trb->t);
1082 }
1083
1084 /* Create, populate and submit the Paired Samples Correlation box */
1085 static void
1086 pscbox (struct t_test_proc *proc)
1087 {
1088   const int rows=1+proc->n_pairs;
1089   const int cols=5;
1090   int i;
1091
1092   struct tab_table *table;
1093
1094   table = tab_create (cols, rows, 0);
1095
1096   tab_columns (table, SOM_COL_DOWN, 1);
1097   tab_headers (table, 0, 0, 1, 0);
1098   tab_box (table, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols - 1, rows - 1);
1099   tab_hline (table, TAL_2, 0, cols - 1, 1);
1100   tab_vline (table, TAL_2, 2, 0, rows - 1);
1101   tab_dim (table, tab_natural_dimensions, NULL);
1102   tab_title (table, _("Paired Samples Correlations"));
1103
1104   /* column headings */
1105   tab_text (table, 2, 0, TAB_CENTER | TAT_TITLE, _("N"));
1106   tab_text (table, 3, 0, TAB_CENTER | TAT_TITLE, _("Correlation"));
1107   tab_text (table, 4, 0, TAB_CENTER | TAT_TITLE, _("Sig."));
1108
1109   for (i = 0; i < proc->n_pairs; i++)
1110     {
1111       struct pair *pair = &proc->pairs[i];
1112       double p, q;
1113       double df = pair->n -2;
1114       double correlation_t = (pair->correlation * sqrt (df) /
1115                               sqrt (1 - pow2 (pair->correlation)));
1116
1117       /* row headings */
1118       tab_text_format (table, 0, i + 1, TAB_LEFT | TAT_TITLE,
1119                        _("Pair %d"), i);
1120       tab_text_format (table, 1, i + 1, TAB_LEFT | TAT_TITLE,
1121                        _("%s & %s"),
1122                        var_get_name (pair->v[0]),
1123                        var_get_name (pair->v[1]));
1124
1125       /* row data */
1126       tab_double (table, 2, i + 1, TAB_RIGHT, pair->n, &proc->weight_format);
1127       tab_double (table, 3, i + 1, TAB_RIGHT, pair->correlation, NULL);
1128
1129       p = gsl_cdf_tdist_P (correlation_t, df);
1130       q = gsl_cdf_tdist_Q (correlation_t, df);
1131       tab_double (table, 4, i + 1, TAB_RIGHT,
1132                  2.0 * (correlation_t > 0 ? q : p), NULL);
1133     }
1134
1135   tab_submit (table);
1136 }
1137 \f
1138 /* Calculation Implementation */
1139
1140 /* Calculations common to all variants of the T test. */
1141 static void
1142 common_calc (const struct dictionary *dict,
1143              struct t_test_proc *proc,
1144              struct casereader *reader)
1145 {
1146   struct ccase *c;
1147   int i;
1148
1149   for (i = 0; i < proc->n_vars; i++)
1150     {
1151       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
1152       gs->sum = 0;
1153       gs->n = 0;
1154       gs->ssq = 0;
1155       gs->sum_diff = 0;
1156     }
1157
1158   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
1159     {
1160       double weight = dict_get_case_weight (dict, c, NULL);
1161
1162       /* Listwise has to be implicit if the independent variable
1163          is missing ?? */
1164       if (proc->mode == T_IND_SAMPLES)
1165         {
1166           if (var_is_value_missing (proc->indep_var,
1167                                     case_data (c, proc->indep_var),
1168                                     proc->exclude))
1169             continue;
1170         }
1171
1172       for (i = 0; i < proc->n_vars; i++)
1173         {
1174           const struct variable *v = proc->vars[i];
1175           const union value *val = case_data (c, v);
1176
1177           if (!var_is_value_missing (v, val, proc->exclude))
1178             {
1179               struct group_statistics *gs;
1180               gs = &group_proc_get (v)->ugs;
1181
1182               gs->n += weight;
1183               gs->sum += weight * val->f;
1184               gs->ssq += weight * pow2 (val->f);
1185             }
1186         }
1187     }
1188   casereader_destroy (reader);
1189
1190   for (i = 0; i < proc->n_vars; i++)
1191     {
1192       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
1193
1194       gs->mean = gs->sum / gs->n;
1195       gs->s_std_dev = sqrt (((gs->ssq / gs->n) - pow2 (gs->mean)));
1196       gs->std_dev = sqrt (gs->n / (gs->n- 1)
1197                           * ((gs->ssq / gs->n) - pow2 (gs->mean)));
1198       gs->se_mean = gs->std_dev / sqrt (gs->n);
1199       gs->mean_diff = gs->sum_diff / gs->n;
1200     }
1201 }
1202
1203 /* Calculations for one sample T test. */
1204 static int
1205 one_sample_calc (const struct dictionary *dict, struct t_test_proc *proc,
1206                  struct casereader *reader)
1207 {
1208   struct ccase *c;
1209   int i;
1210
1211   for (i = 0; i < proc->n_vars; i++)
1212     {
1213       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
1214       gs->sum_diff = 0;
1215     }
1216
1217   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
1218     {
1219       double weight = dict_get_case_weight (dict, c, NULL);
1220       for (i = 0; i < proc->n_vars; i++)
1221         {
1222           const struct variable *v = proc->vars[i];
1223           struct group_statistics *gs = &group_proc_get (v)->ugs;
1224           const union value *val = case_data (c, v);
1225           if (!var_is_value_missing (v, val, proc->exclude))
1226             gs->sum_diff += weight * (val->f - proc->testval);
1227         }
1228     }
1229
1230   for (i = 0; i < proc->n_vars; i++)
1231     {
1232       struct group_statistics *gs = &group_proc_get (proc->vars[i])->ugs;
1233       gs->mean_diff = gs->sum_diff / gs->n;
1234     }
1235
1236   casereader_destroy (reader);
1237
1238   return 0;
1239 }
1240
1241 static int
1242 paired_calc (const struct dictionary *dict, struct t_test_proc *proc,
1243              struct casereader *reader)
1244 {
1245   struct ccase *c;
1246   int i;
1247
1248   for (i = 0; i < proc->n_pairs; i++)
1249     {
1250       struct pair *pair = &proc->pairs[i];
1251       pair->n = 0;
1252       pair->sum[0] = pair->sum[1] = 0;
1253       pair->ssq[0] = pair->ssq[1] = 0;
1254       pair->sum_of_prod = 0;
1255       pair->correlation = 0;
1256       pair->sum_of_diffs = 0;
1257       pair->ssq_diffs = 0;
1258     }
1259
1260   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
1261     {
1262       double weight = dict_get_case_weight (dict, c, NULL);
1263       for (i = 0; i < proc->n_pairs; i++)
1264         {
1265           struct pair *pair = &proc->pairs[i];
1266           const struct variable *v0 = pair->v[0];
1267           const struct variable *v1 = pair->v[1];
1268
1269           const union value *val0 = case_data (c, v0);
1270           const union value *val1 = case_data (c, v1);
1271
1272           if (!var_is_value_missing (v0, val0, proc->exclude)
1273               && !var_is_value_missing (v1, val1, proc->exclude))
1274             {
1275               pair->n += weight;
1276               pair->sum[0] += weight * val0->f;
1277               pair->sum[1] += weight * val1->f;
1278               pair->ssq[0] += weight * pow2 (val0->f);
1279               pair->ssq[1] += weight * pow2 (val1->f);
1280               pair->sum_of_prod += weight * val0->f * val1->f;
1281               pair->sum_of_diffs += weight * (val0->f - val1->f);
1282               pair->ssq_diffs += weight * pow2 (val0->f - val1->f);
1283             }
1284         }
1285     }
1286
1287   for (i = 0; i < proc->n_pairs; i++)
1288     {
1289       struct pair *pair = &proc->pairs[i];
1290       const double n = pair->n;
1291       int j;
1292
1293       for (j=0; j < 2; j++)
1294         {
1295           pair->mean[j] = pair->sum[j] / n;
1296           pair->s_std_dev[j] = sqrt ((pair->ssq[j] / n
1297                                       - pow2 (pair->mean[j])));
1298           pair->std_dev[j] = sqrt (n / (n- 1) * (pair->ssq[j] / n
1299                                                 - pow2 (pair->mean[j])));
1300         }
1301
1302       pair->correlation = (pair->sum_of_prod / pair->n
1303                            - pair->mean[0] * pair->mean[1]);
1304       /* correlation now actually contains the covariance */
1305       pair->correlation /= pair->std_dev[0] * pair->std_dev[1];
1306       pair->correlation *= pair->n / (pair->n - 1);
1307
1308       pair->mean_diff = pair->sum_of_diffs / n;
1309       pair->std_dev_diff = sqrt (n / (n - 1) * ((pair->ssq_diffs / n)
1310                                                 - pow2 (pair->mean_diff)));
1311     }
1312
1313   casereader_destroy (reader);
1314   return 0;
1315 }
1316
1317 static int
1318 group_calc (const struct dictionary *dict, struct t_test_proc *proc,
1319             struct casereader *reader)
1320 {
1321   struct ccase *c;
1322   int i;
1323
1324   for (i = 0; i < proc->n_vars; i++)
1325     {
1326       struct group_proc *ttpr = group_proc_get (proc->vars[i]);
1327       int j;
1328
1329       /* There's always 2 groups for a T - TEST */
1330       ttpr->n_groups = 2;
1331       ttpr->group_hash = hsh_create (2,
1332                                      (hsh_compare_func *) compare_group_binary,
1333                                      (hsh_hash_func *) hash_group_binary,
1334                                      (hsh_free_func *) free_group,
1335                                      proc);
1336
1337       for (j = 0; j < 2; j++)
1338         {
1339           struct group_statistics *gs = xmalloc (sizeof *gs);
1340           gs->sum = 0;
1341           gs->n = 0;
1342           gs->ssq = 0;
1343           if (proc->criterion == CMP_EQ)
1344             gs->id = proc->g_value[j];
1345           else
1346             {
1347               if (j == 0)
1348                 gs->id.f = proc->critical_value - 1.0;
1349               else
1350                 gs->id.f = proc->critical_value + 1.0;
1351             }
1352
1353           hsh_insert (ttpr->group_hash, gs);
1354         }
1355     }
1356
1357   for (; (c = casereader_read (reader)) != NULL; case_unref (c))
1358     {
1359       const double weight = dict_get_case_weight (dict, c, NULL);
1360       const union value *gv;
1361
1362       if (var_is_value_missing (proc->indep_var,
1363                                 case_data (c, proc->indep_var), proc->exclude))
1364         continue;
1365
1366       gv = case_data (c, proc->indep_var);
1367       for (i = 0; i < proc->n_vars; i++)
1368         {
1369           const struct variable *var = proc->vars[i];
1370           const union value *val = case_data (c, var);
1371           struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1372           struct group_statistics *gs = hsh_find (grp_hash, gv);
1373
1374           /* If the independent variable doesn't match either of the values
1375              for this case then move on to the next case. */
1376           if (gs == NULL)
1377             break;
1378
1379           if (!var_is_value_missing (var, val, proc->exclude))
1380             {
1381               gs->n += weight;
1382               gs->sum += weight * val->f;
1383               gs->ssq += weight * pow2 (val->f);
1384             }
1385         }
1386     }
1387
1388   for (i = 0; i < proc->n_vars; i++)
1389     {
1390       const struct variable *var = proc->vars[i];
1391       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1392       struct hsh_iterator g;
1393       struct group_statistics *gs;
1394       int count = 0;
1395
1396       for (gs = hsh_first (grp_hash, &g); gs != NULL;
1397            gs = hsh_next (grp_hash, &g))
1398         {
1399           gs->mean = gs->sum / gs->n;
1400           gs->s_std_dev = sqrt (((gs->ssq / gs->n) - pow2 (gs->mean)));
1401           gs->std_dev = sqrt (gs->n / (gs->n- 1)
1402                               * ((gs->ssq / gs->n) - pow2 (gs->mean)));
1403           gs->se_mean = gs->std_dev / sqrt (gs->n);
1404           count++;
1405         }
1406       assert (count == 2);
1407     }
1408
1409   casereader_destroy (reader);
1410
1411   return 0;
1412 }
1413
1414 static void
1415 calculate (struct t_test_proc *proc,
1416            struct casereader *input, const struct dataset *ds)
1417 {
1418   const struct dictionary *dict = dataset_dict (ds);
1419   struct ssbox stat_summary_box;
1420   struct trbox test_results_box;
1421   struct taint *taint;
1422   struct ccase *c;
1423
1424   c = casereader_peek (input, 0);
1425   if (c == NULL)
1426     {
1427       casereader_destroy (input);
1428       return;
1429     }
1430   output_split_file_values (ds, c);
1431   case_unref (c);
1432
1433   if (proc->listwise_missing)
1434     input = casereader_create_filter_missing (input,
1435                                               proc->vars,
1436                                               proc->n_vars,
1437                                               proc->exclude, NULL, NULL);
1438   input = casereader_create_filter_weight (input, dict, NULL, NULL);
1439   taint = taint_clone (casereader_get_taint (input));
1440
1441   common_calc (dict, proc, casereader_clone (input));
1442   switch (proc->mode)
1443     {
1444     case T_1_SAMPLE:
1445       one_sample_calc (dict, proc, input);
1446       break;
1447     case T_PAIRED:
1448       paired_calc (dict, proc, input);
1449       break;
1450     case T_IND_SAMPLES:
1451       group_calc (dict, proc, casereader_clone (input));
1452       levene (dict, input, proc->indep_var, proc->n_vars, proc->vars,
1453               proc->exclude);
1454       break;
1455     default:
1456       NOT_REACHED ();
1457     }
1458
1459   if (!taint_has_tainted_successor (taint))
1460     {
1461       ssbox_create (&stat_summary_box, proc);
1462       ssbox_populate (&stat_summary_box, proc);
1463       ssbox_finalize (&stat_summary_box);
1464
1465       if (proc->mode == T_PAIRED)
1466         pscbox (proc);
1467
1468       trbox_create (&test_results_box, proc);
1469       trbox_populate (&test_results_box, proc);
1470       trbox_finalize (&test_results_box);
1471     }
1472
1473   taint_destroy (taint);
1474 }
1475
1476 /* return 0 if G belongs to group 0,
1477           1 if it belongs to group 1,
1478           2 if it belongs to neither group */
1479 static int
1480 which_group (const struct group_statistics *g,
1481              const struct t_test_proc *proc)
1482 {
1483   int width = var_get_width (proc->indep_var);
1484
1485   if (0 == value_compare_3way (&g->id, &proc->g_value[0], width))
1486     return 0;
1487
1488   if (0 == value_compare_3way (&g->id, &proc->g_value[1], width))
1489     return 1;
1490
1491   return 2;
1492 }
1493
1494 /* Return -1 if the id of a is less than b; +1 if greater than and
1495    0 if equal */
1496 static int
1497 compare_group_binary (const struct group_statistics *a,
1498                       const struct group_statistics *b,
1499                       const struct t_test_proc *proc)
1500 {
1501   int flag_a;
1502   int flag_b;
1503
1504   if (proc->criterion == CMP_LE)
1505     {
1506       flag_a = (a->id.f < proc->critical_value);
1507       flag_b = (b->id.f < proc->critical_value);
1508     }
1509   else
1510     {
1511       flag_a = which_group (a, proc);
1512       flag_b = which_group (b, proc);
1513     }
1514
1515   if (flag_a < flag_b)
1516     return - 1;
1517
1518   return (flag_a > flag_b);
1519 }
1520
1521 /* This is a degenerate case of a hash, since it can only return three possible
1522    values.  It's really a comparison, being used as a hash function */
1523
1524 static unsigned
1525 hash_group_binary (const struct group_statistics *g,
1526                    const struct t_test_proc *proc)
1527 {
1528   return (proc->criterion == CMP_LE
1529           ? g->id.f < proc->critical_value
1530           : which_group (g, proc));
1531 }
1532
1533 /*
1534   Local Variables:
1535   mode: c
1536   End:
1537 */