pspp: Add --no-output option to allow entirely disabling output.
[pspp] / src / language / stats / t-test.q
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 1997-9, 2000, 2009 Free Software Foundation, Inc.
3
4    This program is free software: you can redistribute it and/or modify
5    it under the terms of the GNU General Public License as published by
6    the Free Software Foundation, either version 3 of the License, or
7    (at your option) any later version.
8
9    This program is distributed in the hope that it will be useful,
10    but WITHOUT ANY WARRANTY; without even the implied warranty of
11    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
12    GNU General Public License for more details.
13
14    You should have received a copy of the GNU General Public License
15    along with this program.  If not, see <http://www.gnu.org/licenses/>. */
16
17 #include <config.h>
18
19 #include <gsl/gsl_cdf.h>
20 #include <math.h>
21 #include <stdint.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24
25 #include <data/case.h>
26 #include <data/casegrouper.h>
27 #include <data/casereader.h>
28 #include <data/dictionary.h>
29 #include <data/procedure.h>
30 #include <data/value-labels.h>
31 #include <data/variable.h>
32 #include <language/command.h>
33 #include <language/dictionary/split-file.h>
34 #include <language/lexer/lexer.h>
35 #include <language/lexer/value-parser.h>
36 #include <libpspp/array.h>
37 #include <libpspp/assertion.h>
38 #include <libpspp/compiler.h>
39 #include <libpspp/hash.h>
40 #include <libpspp/message.h>
41 #include <libpspp/misc.h>
42 #include <libpspp/str.h>
43 #include <libpspp/taint.h>
44 #include <math/group-proc.h>
45 #include <math/levene.h>
46 #include <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, '=');
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, '('))
291     n_values = 0;
292   else
293     {
294       if (!parse_value (lexer, &proc->g_value[0], width))
295         return 0;
296       lex_match (lexer, ',');
297       if (lex_match (lexer, ')'))
298         n_values = 1;
299       else
300         {
301           if (!parse_value (lexer, &proc->g_value[1], width)
302               || !lex_force_match (lexer, ')'))
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, '=');
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, '(')
374           && lex_match_id (lexer, "PAIRED")
375           && lex_match (lexer, ')'))
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 static void
1405 calculate (struct t_test_proc *proc,
1406            struct casereader *input, const struct dataset *ds)
1407 {
1408   const struct dictionary *dict = dataset_dict (ds);
1409   struct ssbox stat_summary_box;
1410   struct trbox test_results_box;
1411   struct taint *taint;
1412   struct ccase *c;
1413
1414   c = casereader_peek (input, 0);
1415   if (c == NULL)
1416     {
1417       casereader_destroy (input);
1418       return;
1419     }
1420   output_split_file_values (ds, c);
1421   case_unref (c);
1422
1423   if (proc->listwise_missing)
1424     input = casereader_create_filter_missing (input,
1425                                               proc->vars,
1426                                               proc->n_vars,
1427                                               proc->exclude, NULL, NULL);
1428   input = casereader_create_filter_weight (input, dict, NULL, NULL);
1429   taint = taint_clone (casereader_get_taint (input));
1430
1431   common_calc (dict, proc, casereader_clone (input));
1432   switch (proc->mode)
1433     {
1434     case T_1_SAMPLE:
1435       one_sample_calc (dict, proc, input);
1436       break;
1437     case T_PAIRED:
1438       paired_calc (dict, proc, input);
1439       break;
1440     case T_IND_SAMPLES:
1441       group_calc (dict, proc, casereader_clone (input));
1442       levene (dict, input, proc->indep_var, proc->n_vars, proc->vars,
1443               proc->exclude);
1444       break;
1445     default:
1446       NOT_REACHED ();
1447     }
1448
1449   if (!taint_has_tainted_successor (taint))
1450     {
1451       ssbox_create (&stat_summary_box, proc);
1452       ssbox_populate (&stat_summary_box, proc);
1453       ssbox_finalize (&stat_summary_box);
1454
1455       if (proc->mode == T_PAIRED)
1456         pscbox (proc);
1457
1458       trbox_create (&test_results_box, proc);
1459       trbox_populate (&test_results_box, proc);
1460       trbox_finalize (&test_results_box);
1461     }
1462
1463   taint_destroy (taint);
1464 }
1465
1466 /* return 0 if G belongs to group 0,
1467           1 if it belongs to group 1,
1468           2 if it belongs to neither group */
1469 static int
1470 which_group (const struct group_statistics *g,
1471              const struct t_test_proc *proc)
1472 {
1473   int width = var_get_width (proc->indep_var);
1474
1475   if (0 == value_compare_3way (&g->id, &proc->g_value[0], width))
1476     return 0;
1477
1478   if (0 == value_compare_3way (&g->id, &proc->g_value[1], width))
1479     return 1;
1480
1481   return 2;
1482 }
1483
1484 /* Return -1 if the id of a is less than b; +1 if greater than and
1485    0 if equal */
1486 static int
1487 compare_group_binary (const struct group_statistics *a,
1488                       const struct group_statistics *b,
1489                       const struct t_test_proc *proc)
1490 {
1491   int flag_a;
1492   int flag_b;
1493
1494   if (proc->criterion == CMP_LE)
1495     {
1496       flag_a = (a->id.f < proc->critical_value);
1497       flag_b = (b->id.f < proc->critical_value);
1498     }
1499   else
1500     {
1501       flag_a = which_group (a, proc);
1502       flag_b = which_group (b, proc);
1503     }
1504
1505   if (flag_a < flag_b)
1506     return - 1;
1507
1508   return (flag_a > flag_b);
1509 }
1510
1511 /* This is a degenerate case of a hash, since it can only return three possible
1512    values.  It's really a comparison, being used as a hash function */
1513
1514 static unsigned
1515 hash_group_binary (const struct group_statistics *g,
1516                    const struct t_test_proc *proc)
1517 {
1518   return (proc->criterion == CMP_LE
1519           ? g->id.f < proc->critical_value
1520           : which_group (g, proc));
1521 }
1522
1523 /*
1524   Local Variables:
1525   mode: c
1526   End:
1527 */