src/language/stats/t-test.q

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