pintos-os.org Git - pspp/blob - src/language/stats/t-test.q

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