8446f8355e50afcf00c8d9826064f2ef317cd498
[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 <libpspp/assertion.h>
36 #include <libpspp/compiler.h>
37 #include <libpspp/hash.h>
38 #include <libpspp/message.h>
39 #include <libpspp/misc.h>
40 #include <libpspp/str.h>
41 #include <libpspp/taint.h>
42 #include <math/group-proc.h>
43 #include <math/levene.h>
44 #include <output/manager.h>
45 #include <output/table.h>
46 #include <data/format.h>
47
48 #include "xalloc.h"
49
50 #include "gettext.h"
51 #define _(msgid) gettext (msgid)
52
53 /* (headers) */
54
55 /* (specification)
56    "T-TEST" (tts_):
57      +groups=custom;
58      testval=double;
59      +variables=varlist("PV_NO_SCRATCH | PV_NUMERIC");
60      +pairs=custom;
61      missing=miss:!analysis/listwise,
62             incl:include/!exclude;
63      +format=fmt:!labels/nolabels;
64      criteria=:cin(d:criteria,"%s > 0. && %s < 1.").
65 */
66 /* (declarations) */
67 /* (functions) */
68
69
70 /* Variable for the GROUPS subcommand, if given. */
71 static struct variable *indep_var;
72
73 enum comparison
74   {
75     CMP_LE = -2,
76     CMP_EQ = 0,
77   };
78
79 struct group_properties
80 {
81   /* The comparison criterion */
82   enum comparison criterion;
83
84   /* The independent variable */
85   struct variable *indep_var;
86
87   union {
88     /* The value of the independent variable at which groups are determined to
89        belong to one group or the other */
90     double critical_value;
91
92
93     /* The values of the independent variable for each group */
94     union value g_value[2];
95   } v ;
96
97 };
98
99
100 static struct group_properties gp ;
101
102
103
104 /* PAIRS: Number of pairs to be compared ; each pair. */
105 static int n_pairs = 0 ;
106 struct pair
107 {
108   /* The variables comprising the pair */
109   const struct variable *v[2];
110
111   /* The number of valid variable pairs */
112   double n;
113
114   /* The sum of the members */
115   double sum[2];
116
117   /* sum of squares of the members */
118   double ssq[2];
119
120   /* Std deviation of the members */
121   double std_dev[2];
122
123
124   /* Sample Std deviation of the members */
125   double s_std_dev[2];
126
127   /* The means of the members */
128   double mean[2];
129
130   /* The correlation coefficient between the variables */
131   double correlation;
132
133   /* The sum of the differences */
134   double sum_of_diffs;
135
136   /* The sum of the products */
137   double sum_of_prod;
138
139   /* The mean of the differences */
140   double mean_diff;
141
142   /* The sum of the squares of the differences */
143   double ssq_diffs;
144
145   /* The std deviation of the differences */
146   double std_dev_diff;
147 };
148
149 static struct pair *pairs=0;
150
151 static int parse_value (struct lexer *lexer, union value * v, enum val_type);
152
153 /* Structures and Functions for the Statistics Summary Box */
154 struct ssbox;
155 typedef void populate_ssbox_func (struct ssbox *ssb,
156                                   const struct dictionary *,
157                                   struct cmd_t_test *cmd);
158 typedef void finalize_ssbox_func (struct ssbox *ssb);
159
160 struct ssbox
161 {
162   struct tab_table *t;
163
164   populate_ssbox_func *populate;
165   finalize_ssbox_func *finalize;
166
167 };
168
169 /* Create a ssbox */
170 void ssbox_create (struct ssbox *ssb,   struct cmd_t_test *cmd, int mode);
171
172 /* Populate a ssbox according to cmd */
173 void ssbox_populate (struct ssbox *ssb, const struct dictionary *dict,
174                      struct cmd_t_test *cmd);
175
176 /* Submit and destroy a ssbox */
177 void ssbox_finalize (struct ssbox *ssb);
178
179 /* A function to create, populate and submit the Paired Samples Correlation
180    box */
181 static void pscbox (const struct dictionary *);
182
183
184 /* Structures and Functions for the Test Results Box */
185 struct trbox;
186
187 typedef void populate_trbox_func (struct trbox *trb,
188                                   const struct dictionary *dict,
189                                   struct cmd_t_test *cmd);
190 typedef void finalize_trbox_func (struct trbox *trb);
191
192 struct trbox {
193   struct tab_table *t;
194   populate_trbox_func *populate;
195   finalize_trbox_func *finalize;
196 };
197
198 /* Create a trbox */
199 void trbox_create (struct trbox *trb,   struct cmd_t_test *cmd, int mode);
200
201 /* Populate a ssbox according to cmd */
202 static void trbox_populate (struct trbox *trb, const struct dictionary *dict,
203                      struct cmd_t_test *cmd);
204
205 /* Submit and destroy a ssbox */
206 void trbox_finalize (struct trbox *trb);
207
208 /* Which mode was T-TEST invoked */
209 enum {
210   T_1_SAMPLE = 0 ,
211   T_IND_SAMPLES,
212   T_PAIRED
213 };
214
215
216 static int common_calc (const struct dictionary *dict,
217                         const struct ccase *, void *,
218                         enum mv_class);
219 static void common_precalc (struct cmd_t_test *);
220 static void common_postcalc (struct cmd_t_test *);
221
222 static int one_sample_calc (const struct dictionary *dict, const struct ccase *, void *, enum mv_class);
223 static void one_sample_precalc (struct cmd_t_test *);
224 static void one_sample_postcalc (struct cmd_t_test *);
225
226 static int  paired_calc (const struct dictionary *dict, const struct ccase *,
227                          struct cmd_t_test*, enum mv_class);
228 static void paired_precalc (struct cmd_t_test *);
229 static void paired_postcalc (struct cmd_t_test *);
230
231 static void group_precalc (struct cmd_t_test *);
232 static int  group_calc (const struct dictionary *dict, const struct ccase *,
233                         struct cmd_t_test *, enum mv_class);
234 static void group_postcalc (struct cmd_t_test *);
235
236
237 static void calculate (struct cmd_t_test *,
238                       struct casereader *,
239                       const struct dataset *);
240
241 static  int mode;
242
243 static int compare_group_binary (const struct group_statistics *a,
244                                 const struct group_statistics *b,
245                                 const struct group_properties *p);
246
247
248 static unsigned  hash_group_binary (const struct group_statistics *g,
249                                    const struct group_properties *p);
250
251
252
253 int
254 cmd_t_test (struct lexer *lexer, struct dataset *ds)
255 {
256   struct cmd_t_test cmd;
257   struct casegrouper *grouper;
258   struct casereader *group;
259   bool ok;
260
261   if ( !parse_t_test (lexer, ds, &cmd, NULL) )
262     return CMD_FAILURE;
263
264   if (! cmd.sbc_criteria)
265     cmd.criteria=0.95;
266
267   {
268     int m=0;
269     if (cmd.sbc_testval) ++m;
270     if (cmd.sbc_groups) ++m;
271     if (cmd.sbc_pairs) ++m;
272
273     if ( m != 1)
274       {
275         msg (SE,
276             _ ("TESTVAL, GROUPS and PAIRS subcommands are mutually exclusive.")
277             );
278         free_t_test (&cmd);
279         return CMD_FAILURE;
280       }
281   }
282
283   if (cmd.sbc_testval)
284     mode=T_1_SAMPLE;
285   else if (cmd.sbc_groups)
286     mode=T_IND_SAMPLES;
287   else
288     mode=T_PAIRED;
289
290   if ( mode == T_PAIRED)
291     {
292       if (cmd.sbc_variables)
293         {
294           msg (SE, _ ("VARIABLES subcommand is not appropriate with PAIRS"));
295           free_t_test (&cmd);
296           return CMD_FAILURE;
297         }
298       else
299         {
300           /* Iterate through the pairs and put each variable that is a
301              member of a pair into cmd.v_variables */
302
303           int i;
304           struct hsh_iterator hi;
305           struct const_hsh_table *hash;
306           const struct variable *v;
307
308           hash = const_hsh_create (n_pairs, compare_vars_by_name, hash_var_by_name,
309           0, 0);
310
311           for (i=0; i < n_pairs; ++i)
312             {
313               const_hsh_insert (hash, pairs[i].v[0]);
314               const_hsh_insert (hash, pairs[i].v[1]);
315             }
316
317           assert (cmd.n_variables == 0);
318           cmd.n_variables = const_hsh_count (hash);
319
320           cmd.v_variables = xnrealloc (cmd.v_variables, cmd.n_variables,
321                                        sizeof *cmd.v_variables);
322           /* Iterate through the hash */
323           for (i=0,v = const_hsh_first (hash, &hi);
324                v != 0;
325                v = const_hsh_next (hash, &hi) )
326             cmd.v_variables[i++]=v;
327           const_hsh_destroy (hash);
328         }
329     }
330   else if ( !cmd.sbc_variables)
331     {
332       msg (SE, _ ("One or more VARIABLES must be specified."));
333       free_t_test (&cmd);
334       return CMD_FAILURE;
335     }
336
337   /* Data pass. */
338   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
339   while (casegrouper_get_next_group (grouper, &group))
340     calculate (&cmd, group, ds);
341   ok = casegrouper_destroy (grouper);
342   ok = proc_commit (ds) && ok;
343
344   n_pairs=0;
345   free (pairs);
346   pairs=0;
347
348   if ( mode == T_IND_SAMPLES)
349     {
350       int v;
351       /* Destroy any group statistics we created */
352       for (v = 0 ; v < cmd.n_variables ; ++v )
353         {
354           struct group_proc *grpp = group_proc_get (cmd.v_variables[v]);
355           hsh_destroy (grpp->group_hash);
356         }
357     }
358
359   free_t_test (&cmd);
360   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
361 }
362
363 static int
364 tts_custom_groups (struct lexer *lexer, struct dataset *ds, struct cmd_t_test *cmd UNUSED, 
365         void *aux UNUSED)
366 {
367   int n_group_values=0;
368
369   lex_match (lexer, '=');
370
371   indep_var = parse_variable (lexer, dataset_dict (ds));
372   if (!indep_var)
373     {
374       lex_error (lexer, "expecting variable name in GROUPS subcommand");
375       return 0;
376     }
377
378   if (var_is_long_string (indep_var))
379     {
380       msg (SE, _ ("Long string variable %s is not valid here."),
381            var_get_name (indep_var));
382       return 0;
383     }
384
385   if (!lex_match (lexer, '('))
386     {
387       if (var_is_numeric (indep_var))
388         {
389           gp.v.g_value[0].f = 1;
390           gp.v.g_value[1].f = 2;
391
392           gp.criterion = CMP_EQ;
393
394           n_group_values = 2;
395
396           return 1;
397         }
398       else
399         {
400           msg (SE, _ ("When applying GROUPS to a string variable, two "
401                      "values must be specified."));
402           return 0;
403         }
404     }
405
406   if (!parse_value (lexer, &gp.v.g_value[0], var_get_width (indep_var)))
407       return 0;
408
409   lex_match (lexer, ',');
410   if (lex_match (lexer, ')'))
411     {
412       if (var_is_alpha (indep_var))
413         {
414           msg (SE, _ ("When applying GROUPS to a string variable, two "
415                      "values must be specified."));
416           return 0;
417         }
418       gp.criterion = CMP_LE;
419       gp.v.critical_value = gp.v.g_value[0].f;
420
421       n_group_values = 1;
422       return 1;
423     }
424
425   if (!parse_value (lexer, &gp.v.g_value[1], var_get_width (indep_var)))
426     return 0;
427
428   n_group_values = 2;
429   if (!lex_force_match (lexer, ')'))
430     return 0;
431
432   if ( n_group_values == 2 )
433     gp.criterion = CMP_EQ ;
434   else
435     gp.criterion = CMP_LE ;
436
437
438   if ( var_is_alpha (indep_var))
439     {
440       buf_copy_rpad (gp.v.g_value [0].s, var_get_width (indep_var),
441                      gp.v.g_value [0].s, strlen (gp.v.g_value[0].s));
442
443       buf_copy_rpad (gp.v.g_value [1].s, var_get_width (indep_var),
444                      gp.v.g_value [1].s, strlen (gp.v.g_value[1].s));
445     }
446
447   return 1;
448 }
449
450
451 static int
452 tts_custom_pairs (struct lexer *lexer, struct dataset *ds, struct cmd_t_test *cmd UNUSED, void *aux UNUSED)
453 {
454   const struct variable **vars;
455   size_t n_vars;
456   size_t n_pairs_local;
457
458   size_t n_before_WITH;
459   size_t n_after_WITH = SIZE_MAX;
460   int paired ; /* Was the PAIRED keyword given ? */
461
462   lex_match (lexer, '=');
463
464   n_vars=0;
465   if (!parse_variables_const (lexer, dataset_dict (ds), &vars, &n_vars,
466                         PV_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH))
467     {
468       free (vars);
469       return 0;
470     }
471   assert (n_vars);
472
473   n_before_WITH = 0;
474   if (lex_match (lexer, T_WITH))
475     {
476       n_before_WITH = n_vars;
477       if (!parse_variables_const (lexer, dataset_dict (ds), &vars, &n_vars,
478                             PV_DUPLICATE | PV_APPEND
479                             | PV_NUMERIC | PV_NO_SCRATCH))
480         {
481           free (vars);
482           return 0;
483         }
484       n_after_WITH = n_vars - n_before_WITH;
485     }
486
487   paired = (lex_match (lexer, '(') && lex_match_id (lexer, "PAIRED") && lex_match (lexer, ')'));
488
489   /* Determine the number of pairs needed */
490   if (paired)
491     {
492       if (n_before_WITH != n_after_WITH)
493         {
494           free (vars);
495           msg (SE, _ ("PAIRED was specified but the number of variables "
496                      "preceding WITH (%zu) did not match the number "
497                      "following (%zu)."),
498                n_before_WITH, n_after_WITH);
499           return 0;
500         }
501       n_pairs_local = n_before_WITH;
502     }
503   else if (n_before_WITH > 0) /* WITH keyword given, but not PAIRED keyword */
504     {
505       n_pairs_local = n_before_WITH * n_after_WITH ;
506     }
507   else /* Neither WITH nor PAIRED keyword given */
508     {
509       if (n_vars < 2)
510         {
511           free (vars);
512           msg (SE, _ ("At least two variables must be specified "
513                      "on PAIRS."));
514           return 0;
515         }
516
517       /* how many ways can you pick 2 from n_vars ? */
518       n_pairs_local = n_vars * (n_vars - 1) / 2;
519     }
520
521
522   /* Allocate storage for the pairs */
523   pairs = xnrealloc (pairs, n_pairs + n_pairs_local, sizeof *pairs);
524
525   /* Populate the pairs with the appropriate variables */
526   if ( paired )
527     {
528       int i;
529
530       assert (n_pairs_local == n_vars / 2);
531       for (i = 0; i < n_pairs_local; ++i)
532         {
533           pairs[i].v[n_pairs] = vars[i];
534           pairs[i].v[n_pairs + 1] = vars[i + n_pairs_local];
535         }
536     }
537   else if (n_before_WITH > 0) /* WITH keyword given, but not PAIRED keyword */
538     {
539       int i,j;
540       size_t p = n_pairs;
541
542       for (i=0 ; i < n_before_WITH ; ++i )
543         {
544           for (j=0 ; j < n_after_WITH ; ++j)
545             {
546               pairs[p].v[0] = vars[i];
547               pairs[p].v[1] = vars[j+n_before_WITH];
548               ++p;
549             }
550         }
551     }
552   else /* Neither WITH nor PAIRED given */
553     {
554       size_t i,j;
555       size_t p=n_pairs;
556
557       for (i=0 ; i < n_vars ; ++i )
558         {
559           for (j=i+1 ; j < n_vars ; ++j)
560             {
561               pairs[p].v[0] = vars[i];
562               pairs[p].v[1] = vars[j];
563               ++p;
564             }
565         }
566     }
567
568   n_pairs+=n_pairs_local;
569
570   free (vars);
571   return 1;
572 }
573
574 /* Parses the current token (numeric or string, depending on type)
575     value v and returns success. */
576 static int
577 parse_value (struct lexer *lexer, union value * v, enum val_type type)
578 {
579   if (type == VAL_NUMERIC)
580     {
581       if (!lex_force_num (lexer))
582         return 0;
583       v->f = lex_tokval (lexer);
584     }
585   else
586     {
587       if (!lex_force_string (lexer))
588         return 0;
589       memset  (v->s, ' ', MAX_SHORT_STRING);
590       strncpy (v->s, ds_cstr (lex_tokstr (lexer)), ds_length (lex_tokstr (lexer)));
591     }
592
593   lex_get (lexer);
594
595   return 1;
596 }
597
598
599 /* Implementation of the SSBOX object */
600
601 void ssbox_base_init (struct ssbox *this, int cols,int rows);
602
603 void ssbox_base_finalize (struct ssbox *ssb);
604
605 void ssbox_one_sample_init (struct ssbox *this,
606                            struct cmd_t_test *cmd );
607
608 void ssbox_independent_samples_init (struct ssbox *this,
609                                     struct cmd_t_test *cmd);
610
611 void ssbox_paired_init (struct ssbox *this,
612                            struct cmd_t_test *cmd);
613
614
615 /* Factory to create an ssbox */
616 void
617 ssbox_create (struct ssbox *ssb, struct cmd_t_test *cmd, int mode)
618 {
619     switch (mode)
620       {
621       case T_1_SAMPLE:
622         ssbox_one_sample_init (ssb,cmd);
623         break;
624       case T_IND_SAMPLES:
625         ssbox_independent_samples_init (ssb,cmd);
626         break;
627       case T_PAIRED:
628         ssbox_paired_init (ssb,cmd);
629         break;
630       default:
631         NOT_REACHED ();
632       }
633 }
634
635
636
637 /* Despatcher for the populate method */
638 void
639 ssbox_populate (struct ssbox *ssb, const struct dictionary *dict,
640                 struct cmd_t_test *cmd)
641 {
642   ssb->populate (ssb, dict, cmd);
643 }
644
645
646 /* Despatcher for finalize */
647 void
648 ssbox_finalize (struct ssbox *ssb)
649 {
650   ssb->finalize (ssb);
651 }
652
653
654 /* Submit the box and clear up */
655 void
656 ssbox_base_finalize (struct ssbox *ssb)
657 {
658   tab_submit (ssb->t);
659 }
660
661
662
663 /* Initialize a ssbox struct */
664 void
665 ssbox_base_init (struct ssbox *this, int cols,int rows)
666 {
667   this->finalize = ssbox_base_finalize;
668   this->t = tab_create (cols, rows, 0);
669
670   tab_columns (this->t, SOM_COL_DOWN, 1);
671   tab_headers (this->t,0,0,1,0);
672   tab_box (this->t, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols -1, rows -1 );
673   tab_hline (this->t, TAL_2,0,cols-1,1);
674   tab_dim (this->t, tab_natural_dimensions);
675 }
676
677 void  ssbox_one_sample_populate (struct ssbox *ssb,
678                                  const struct dictionary *,
679                                  struct cmd_t_test *cmd);
680
681 /* Initialize the one_sample ssbox */
682 void
683 ssbox_one_sample_init (struct ssbox *this,
684                            struct cmd_t_test *cmd )
685 {
686   const int hsize=5;
687   const int vsize=cmd->n_variables+1;
688
689   this->populate = ssbox_one_sample_populate;
690
691   ssbox_base_init (this, hsize,vsize);
692   tab_title (this->t, _ ("One-Sample Statistics"));
693   tab_vline (this->t, TAL_2, 1,0,vsize - 1);
694   tab_text (this->t, 1, 0, TAB_CENTER | TAT_TITLE, _ ("N"));
695   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _ ("Mean"));
696   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _ ("Std. Deviation"));
697   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _ ("SE. Mean"));
698 }
699
700 static void ssbox_independent_samples_populate (struct ssbox *ssb,
701                                                 const struct dictionary *,
702                                                 struct cmd_t_test *cmd);
703
704 /* Initialize the independent samples ssbox */
705 void
706 ssbox_independent_samples_init (struct ssbox *this,
707         struct cmd_t_test *cmd)
708 {
709   int hsize=6;
710   int vsize = cmd->n_variables*2 +1;
711
712   this->populate = ssbox_independent_samples_populate;
713
714   ssbox_base_init (this, hsize,vsize);
715   tab_vline (this->t, TAL_GAP, 1, 0,vsize - 1);
716   tab_title (this->t, _ ("Group Statistics"));
717   tab_text (this->t, 1, 0, TAB_CENTER | TAT_TITLE, var_get_name (indep_var));
718   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _ ("N"));
719   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _ ("Mean"));
720   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _ ("Std. Deviation"));
721   tab_text (this->t, 5, 0, TAB_CENTER | TAT_TITLE, _ ("SE. Mean"));
722 }
723
724
725 /* Populate the ssbox for independent samples */
726 static void
727 ssbox_independent_samples_populate (struct ssbox *ssb,
728                                     const struct dictionary *dict,
729                                     struct cmd_t_test *cmd)
730 {
731   int i;
732
733   const struct variable *wv = dict_get_weight (dict);
734   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
735
736   char *val_lab[2] = {NULL, NULL};
737   double indep_value[2];
738
739   char prefix[2][3]={"",""};
740
741   if ( var_is_numeric (indep_var) )
742     {
743       const char *s;
744
745       s = var_lookup_value_label (indep_var, &gp.v.g_value[0]);
746       val_lab[0] = s ? xstrdup (s) : NULL;
747
748       s = var_lookup_value_label (indep_var, &gp.v.g_value[1]);
749       val_lab[1] = s ? xstrdup (s) : NULL;
750     }
751   else
752     {
753       val_lab[0] = xcalloc (sizeof (char), MAX_SHORT_STRING + 1);
754       val_lab[1] = xcalloc (sizeof (char), MAX_SHORT_STRING + 1);
755       memcpy (val_lab[0], gp.v.g_value[0].s, MAX_SHORT_STRING);
756       memcpy (val_lab[1], gp.v.g_value[1].s, MAX_SHORT_STRING);
757     }
758
759   if (gp.criterion == CMP_LE )
760     {
761       strcpy (prefix[0],">=");
762       strcpy (prefix[1],"<");
763       indep_value[0] = gp.v.critical_value;
764       indep_value[1] = gp.v.critical_value;
765     }
766   else
767     {
768       indep_value[0] = gp.v.g_value[0].f;
769       indep_value[1] = gp.v.g_value[1].f;
770     }
771
772   assert (ssb->t);
773
774   for (i=0; i < cmd->n_variables; ++i)
775     {
776       const struct variable *var = cmd->v_variables[i];
777       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
778       int count=0;
779
780       tab_text (ssb->t, 0, i*2+1, TAB_LEFT,
781                 var_get_name (cmd->v_variables[i]));
782
783       if (val_lab[0])
784         tab_text (ssb->t, 1, i*2+1, TAB_LEFT | TAT_PRINTF,
785                   "%s%s", prefix[0], val_lab[0]);
786       else
787           tab_text (ssb->t, 1, i*2+1, TAB_LEFT | TAT_PRINTF,
788                     "%s%g", prefix[0], indep_value[0]);
789
790
791       if (val_lab[1])
792         tab_text (ssb->t, 1, i*2+1+1, TAB_LEFT | TAT_PRINTF,
793                   "%s%s", prefix[1], val_lab[1]);
794       else
795           tab_text (ssb->t, 1, i*2+1+1, TAB_LEFT | TAT_PRINTF,
796                     "%s%g", prefix[1], indep_value[1]);
797
798
799       /* Fill in the group statistics */
800       for ( count = 0 ; count < 2 ; ++count )
801         {
802           union value search_val;
803
804           struct group_statistics *gs;
805
806           if ( gp.criterion == CMP_LE )
807             {
808               if ( count == 0 )
809                 {
810                   /* >= case  */
811                   search_val.f = gp.v.critical_value + 1.0;
812                 }
813               else
814                 {
815                   /*  less than ( < )  case */
816                   search_val.f = gp.v.critical_value - 1.0;
817                 }
818             }
819           else
820             {
821               search_val = gp.v.g_value[count];
822             }
823
824           gs = hsh_find (grp_hash, (void *) &search_val);
825           assert (gs);
826
827           tab_double (ssb->t, 2, i*2+count+1, TAB_RIGHT, gs->n, wfmt);
828           tab_double (ssb->t, 3, i*2+count+1, TAB_RIGHT, gs->mean, NULL);
829           tab_double (ssb->t, 4, i*2+count+1, TAB_RIGHT, gs->std_dev, NULL);
830           tab_double (ssb->t, 5, i*2+count+1, TAB_RIGHT, gs->se_mean, NULL);
831         }
832     }
833   free (val_lab[0]);
834   free (val_lab[1]);
835 }
836
837
838 static void ssbox_paired_populate (struct ssbox *ssb,
839                                    const struct dictionary *dict,
840                                    struct cmd_t_test *cmd);
841
842 /* Initialize the paired values ssbox */
843 void
844 ssbox_paired_init (struct ssbox *this, struct cmd_t_test *cmd UNUSED)
845 {
846   int hsize=6;
847
848   int vsize = n_pairs*2+1;
849
850   this->populate = ssbox_paired_populate;
851
852   ssbox_base_init (this, hsize,vsize);
853   tab_title (this->t, _ ("Paired Sample Statistics"));
854   tab_vline (this->t,TAL_GAP,1,0,vsize-1);
855   tab_vline (this->t,TAL_2,2,0,vsize-1);
856   tab_text (this->t, 2, 0, TAB_CENTER | TAT_TITLE, _ ("Mean"));
857   tab_text (this->t, 3, 0, TAB_CENTER | TAT_TITLE, _ ("N"));
858   tab_text (this->t, 4, 0, TAB_CENTER | TAT_TITLE, _ ("Std. Deviation"));
859   tab_text (this->t, 5, 0, TAB_CENTER | TAT_TITLE, _ ("SE. Mean"));
860 }
861
862
863 /* Populate the ssbox for paired values */
864 void
865 ssbox_paired_populate (struct ssbox *ssb, const struct dictionary *dict,
866                        struct cmd_t_test *cmd UNUSED)
867 {
868   int i;
869
870   const struct variable *wv = dict_get_weight (dict);
871   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
872
873   assert (ssb->t);
874
875   for (i=0; i < n_pairs; ++i)
876     {
877       int j;
878
879       tab_text (ssb->t, 0, i*2+1, TAB_LEFT | TAT_PRINTF , _ ("Pair %d"),i);
880
881       for (j=0 ; j < 2 ; ++j)
882         {
883           struct group_statistics *gs;
884
885           gs = &group_proc_get (pairs[i].v[j])->ugs;
886
887           /* Titles */
888
889           tab_text (ssb->t, 1, i*2+j+1, TAB_LEFT,
890                     var_get_name (pairs[i].v[j]));
891
892           /* Values */
893           tab_double (ssb->t,2, i*2+j+1, TAB_RIGHT, pairs[i].mean[j], NULL);
894           tab_double (ssb->t,3, i*2+j+1, TAB_RIGHT, pairs[i].n, wfmt);
895           tab_double (ssb->t,4, i*2+j+1, TAB_RIGHT, pairs[i].std_dev[j], NULL);
896           tab_double (ssb->t,5, i*2+j+1, TAB_RIGHT,
897                       pairs[i].std_dev[j]/sqrt (pairs[i].n), NULL);
898
899         }
900     }
901 }
902
903 /* Populate the one sample ssbox */
904 void
905 ssbox_one_sample_populate (struct ssbox *ssb, const struct dictionary *dict,
906                            struct cmd_t_test *cmd)
907 {
908   int i;
909
910   const struct variable *wv = dict_get_weight (dict);
911   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
912
913   assert (ssb->t);
914
915   for (i=0; i < cmd->n_variables; ++i)
916     {
917       struct group_statistics *gs = &group_proc_get (cmd->v_variables[i])->ugs;
918
919       tab_text (ssb->t, 0, i+1, TAB_LEFT, var_get_name (cmd->v_variables[i]));
920       tab_double (ssb->t,1, i+1, TAB_RIGHT, gs->n, wfmt);
921       tab_double (ssb->t,2, i+1, TAB_RIGHT, gs->mean, NULL);
922       tab_double (ssb->t,3, i+1, TAB_RIGHT, gs->std_dev, NULL);
923       tab_double (ssb->t,4, i+1, TAB_RIGHT, gs->se_mean, NULL);
924     }
925 }
926
927
928
929 /* Implementation of the Test Results box struct */
930
931 void trbox_base_init (struct trbox *self,size_t n_vars, int cols);
932 void trbox_base_finalize (struct trbox *trb);
933
934 void trbox_independent_samples_init (struct trbox *trb,
935                                     struct cmd_t_test *cmd );
936
937 static void trbox_independent_samples_populate (struct trbox *trb,
938                                          const struct dictionary *dict,
939                                          struct cmd_t_test *cmd);
940
941 void trbox_one_sample_init (struct trbox *self,
942                       struct cmd_t_test *cmd );
943
944 static void trbox_one_sample_populate (struct trbox *trb,
945                                 const struct dictionary *,
946                                 struct cmd_t_test *cmd);
947
948 void trbox_paired_init (struct trbox *self,
949                        struct cmd_t_test *cmd );
950
951 static void trbox_paired_populate (struct trbox *trb,
952                                    const struct dictionary *,
953                                    struct cmd_t_test *cmd);
954
955
956
957 /* Create a trbox according to mode*/
958 void
959 trbox_create (struct trbox *trb,
960              struct cmd_t_test *cmd, int mode)
961 {
962     switch (mode)
963       {
964       case T_1_SAMPLE:
965         trbox_one_sample_init (trb,cmd);
966         break;
967       case T_IND_SAMPLES:
968         trbox_independent_samples_init (trb,cmd);
969         break;
970       case T_PAIRED:
971         trbox_paired_init (trb,cmd);
972         break;
973       default:
974         NOT_REACHED ();
975       }
976 }
977
978 /* Populate a trbox according to cmd */
979 static void
980 trbox_populate (struct trbox *trb, const struct dictionary *dict,
981                 struct cmd_t_test *cmd)
982 {
983   trb->populate (trb, dict, cmd);
984 }
985
986 /* Submit and destroy a trbox */
987 void
988 trbox_finalize (struct trbox *trb)
989 {
990   trb->finalize (trb);
991 }
992
993 /* Initialize the independent samples trbox */
994 void
995 trbox_independent_samples_init (struct trbox *self,
996                            struct cmd_t_test *cmd UNUSED)
997 {
998   const int hsize=11;
999   const int vsize=cmd->n_variables*2+3;
1000
1001   assert (self);
1002   self->populate = trbox_independent_samples_populate;
1003
1004   trbox_base_init (self,cmd->n_variables*2,hsize);
1005   tab_title (self->t,_ ("Independent Samples Test"));
1006   tab_hline (self->t,TAL_1,2,hsize-1,1);
1007   tab_vline (self->t,TAL_2,2,0,vsize-1);
1008   tab_vline (self->t,TAL_1,4,0,vsize-1);
1009   tab_box (self->t,-1,-1,-1,TAL_1, 2,1,hsize-2,vsize-1);
1010   tab_hline (self->t,TAL_1, hsize-2,hsize-1,2);
1011   tab_box (self->t,-1,-1,-1,TAL_1, hsize-2,2,hsize-1,vsize-1);
1012   tab_joint_text (self->t, 2, 0, 3, 0,
1013                  TAB_CENTER,_ ("Levene's Test for Equality of Variances"));
1014   tab_joint_text (self->t, 4,0,hsize-1,0,
1015                  TAB_CENTER,_ ("t-test for Equality of Means"));
1016
1017   tab_text (self->t,2,2, TAB_CENTER | TAT_TITLE,_ ("F"));
1018   tab_text (self->t,3,2, TAB_CENTER | TAT_TITLE,_ ("Sig."));
1019   tab_text (self->t,4,2, TAB_CENTER | TAT_TITLE,_ ("t"));
1020   tab_text (self->t,5,2, TAB_CENTER | TAT_TITLE,_ ("df"));
1021   tab_text (self->t,6,2, TAB_CENTER | TAT_TITLE,_ ("Sig. (2-tailed)"));
1022   tab_text (self->t,7,2, TAB_CENTER | TAT_TITLE,_ ("Mean Difference"));
1023   tab_text (self->t,8,2, TAB_CENTER | TAT_TITLE,_ ("Std. Error Difference"));
1024   tab_text (self->t,9,2, TAB_CENTER | TAT_TITLE,_ ("Lower"));
1025   tab_text (self->t,10,2, TAB_CENTER | TAT_TITLE,_ ("Upper"));
1026
1027   tab_joint_text (self->t, 9, 1, 10, 1, TAB_CENTER | TAT_PRINTF,
1028                  _ ("%g%% Confidence Interval of the Difference"),
1029                  cmd->criteria*100.0);
1030
1031 }
1032
1033 /* Populate the independent samples trbox */
1034 static void
1035 trbox_independent_samples_populate (struct trbox *self,
1036                                     const struct dictionary *dict UNUSED,
1037                                     struct cmd_t_test *cmd)
1038 {
1039   int i;
1040
1041   assert (self);
1042   for (i=0; i < cmd->n_variables; ++i)
1043     {
1044       double p,q;
1045
1046       double t;
1047       double df;
1048
1049       double df1, df2;
1050
1051       double pooled_variance;
1052       double std_err_diff;
1053       double mean_diff;
1054
1055       const struct variable *var = cmd->v_variables[i];
1056       struct group_proc *grp_data = group_proc_get (var);
1057
1058       struct hsh_table *grp_hash = grp_data->group_hash;
1059
1060       struct group_statistics *gs0 ;
1061       struct group_statistics *gs1 ;
1062
1063       union value search_val;
1064
1065       if ( gp.criterion == CMP_LE )
1066         search_val.f = gp.v.critical_value - 1.0;
1067       else
1068         search_val = gp.v.g_value[0];
1069
1070       gs0 = hsh_find (grp_hash, (void *) &search_val);
1071       assert (gs0);
1072
1073       if ( gp.criterion == CMP_LE )
1074         search_val.f = gp.v.critical_value + 1.0;
1075       else
1076         search_val = gp.v.g_value[1];
1077
1078       gs1 = hsh_find (grp_hash, (void *) &search_val);
1079       assert (gs1);
1080
1081
1082       tab_text (self->t, 0, i*2+3, TAB_LEFT, var_get_name (cmd->v_variables[i]));
1083
1084       tab_text (self->t, 1, i*2+3, TAB_LEFT, _ ("Equal variances assumed"));
1085
1086
1087       tab_double (self->t, 2, i*2+3, TAB_CENTER, grp_data->levene, NULL);
1088
1089       /* Now work out the significance of the Levene test */
1090       df1 = 1; df2 = grp_data->ugs.n - 2;
1091       q = gsl_cdf_fdist_Q (grp_data->levene, df1, df2);
1092
1093       tab_double (self->t, 3, i*2+3, TAB_CENTER, q, NULL);
1094
1095       df = gs0->n + gs1->n - 2.0 ;
1096       tab_double (self->t, 5, i*2+3, TAB_RIGHT, df, NULL);
1097
1098       pooled_variance = ( (gs0->n )*pow2 (gs0->s_std_dev)
1099                           +
1100                           (gs1->n )*pow2 (gs1->s_std_dev)
1101                         ) / df  ;
1102
1103       t = (gs0->mean - gs1->mean) / sqrt (pooled_variance) ;
1104       t /= sqrt ((gs0->n + gs1->n)/ (gs0->n*gs1->n));
1105
1106       tab_double (self->t, 4, i*2+3, TAB_RIGHT, t, NULL);
1107
1108       p = gsl_cdf_tdist_P (t, df);
1109       q = gsl_cdf_tdist_Q (t, df);
1110
1111       tab_double (self->t, 6, i*2+3, TAB_RIGHT, 2.0* (t>0?q:p), NULL);
1112
1113       mean_diff = gs0->mean - gs1->mean;
1114       tab_double (self->t, 7, i*2+3, TAB_RIGHT, mean_diff, NULL);
1115
1116
1117       std_err_diff = sqrt ( pow2 (gs0->se_mean) + pow2 (gs1->se_mean));
1118       tab_double (self->t, 8, i*2+3, TAB_RIGHT, std_err_diff, NULL);
1119
1120
1121       /* Now work out the confidence interval */
1122       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1123
1124       t = gsl_cdf_tdist_Qinv (q,df);
1125       tab_double (self->t, 9, i*2+3, TAB_RIGHT,
1126                 mean_diff - t * std_err_diff, NULL);
1127
1128       tab_double (self->t, 10, i*2+3, TAB_RIGHT,
1129                 mean_diff + t * std_err_diff, NULL);
1130
1131
1132       {
1133         double se2;
1134       /* Now for the \sigma_1 != \sigma_2 case */
1135       tab_text (self->t, 1, i*2+3+1,
1136                 TAB_LEFT, _ ("Equal variances not assumed"));
1137
1138
1139       se2 = (pow2 (gs0->s_std_dev)/ (gs0->n -1) ) +
1140          (pow2 (gs1->s_std_dev)/ (gs1->n -1) );
1141
1142       t = mean_diff / sqrt (se2) ;
1143       tab_double (self->t, 4, i*2+3+1, TAB_RIGHT, t, NULL);
1144
1145       df = pow2 (se2) / (
1146                        (pow2 (pow2 (gs0->s_std_dev)/ (gs0->n - 1 ))
1147                         / (gs0->n -1 )
1148                         )
1149                        +
1150                        (pow2 (pow2 (gs1->s_std_dev)/ (gs1->n - 1 ))
1151                         / (gs1->n -1 )
1152                         )
1153                        ) ;
1154
1155       tab_double (self->t, 5, i*2+3+1, TAB_RIGHT, df, NULL);
1156
1157       p = gsl_cdf_tdist_P (t, df);
1158       q = gsl_cdf_tdist_Q (t, df);
1159
1160       tab_double (self->t, 6, i*2+3+1, TAB_RIGHT, 2.0* (t>0?q:p), NULL);
1161
1162       /* Now work out the confidence interval */
1163       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1164
1165       t = gsl_cdf_tdist_Qinv (q, df);
1166
1167       tab_double (self->t, 7, i*2+3+1, TAB_RIGHT, mean_diff, NULL);
1168
1169
1170       tab_double (self->t, 8, i*2+3+1, TAB_RIGHT, std_err_diff, NULL);
1171
1172
1173       tab_double (self->t, 9, i*2+3+1, TAB_RIGHT,
1174                 mean_diff - t * std_err_diff, NULL);
1175
1176       tab_double (self->t, 10, i*2+3+1, TAB_RIGHT,
1177                 mean_diff + t * std_err_diff, NULL);
1178       }
1179     }
1180 }
1181
1182 /* Initialize the paired samples trbox */
1183 void
1184 trbox_paired_init (struct trbox *self,
1185                            struct cmd_t_test *cmd UNUSED)
1186 {
1187
1188   const int hsize=10;
1189   const int vsize=n_pairs+3;
1190
1191   self->populate = trbox_paired_populate;
1192
1193   trbox_base_init (self,n_pairs,hsize);
1194   tab_title (self->t, _ ("Paired Samples Test"));
1195   tab_hline (self->t,TAL_1,2,6,1);
1196   tab_vline (self->t,TAL_2,2,0,vsize - 1);
1197   tab_joint_text (self->t,2,0,6,0,TAB_CENTER,_ ("Paired Differences"));
1198   tab_box (self->t,-1,-1,-1,TAL_1, 2,1,6,vsize-1);
1199   tab_box (self->t,-1,-1,-1,TAL_1, 6,0,hsize-1,vsize-1);
1200   tab_hline (self->t,TAL_1,5,6, 2);
1201   tab_vline (self->t,TAL_GAP,6,0,1);
1202
1203   tab_joint_text (self->t, 5, 1, 6, 1, TAB_CENTER | TAT_PRINTF,
1204                  _ ("%g%% Confidence Interval of the Difference"),
1205                  cmd->criteria*100.0);
1206
1207   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _ ("Mean"));
1208   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _ ("Std. Deviation"));
1209   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _ ("Std. Error Mean"));
1210   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _ ("Lower"));
1211   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _ ("Upper"));
1212   tab_text (self->t, 7, 2, TAB_CENTER | TAT_TITLE, _ ("t"));
1213   tab_text (self->t, 8, 2, TAB_CENTER | TAT_TITLE, _ ("df"));
1214   tab_text (self->t, 9, 2, TAB_CENTER | TAT_TITLE, _ ("Sig. (2-tailed)"));
1215 }
1216
1217 /* Populate the paired samples trbox */
1218 static void
1219 trbox_paired_populate (struct trbox *trb,
1220                        const struct dictionary *dict,
1221                        struct cmd_t_test *cmd UNUSED)
1222 {
1223   int i;
1224
1225   const struct variable *wv = dict_get_weight (dict);
1226   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
1227
1228   for (i=0; i < n_pairs; ++i)
1229     {
1230       double p,q;
1231       double se_mean;
1232
1233       double n = pairs[i].n;
1234       double t;
1235       double df = n - 1;
1236
1237       tab_text (trb->t, 0, i+3, TAB_LEFT | TAT_PRINTF, _ ("Pair %d"),i);
1238
1239       tab_text (trb->t, 1, i+3, TAB_LEFT | TAT_PRINTF, "%s - %s",
1240                 var_get_name (pairs[i].v[0]),
1241                 var_get_name (pairs[i].v[1]));
1242
1243       tab_double (trb->t, 2, i+3, TAB_RIGHT, pairs[i].mean_diff, NULL);
1244
1245       tab_double (trb->t, 3, i+3, TAB_RIGHT, pairs[i].std_dev_diff, NULL);
1246
1247       /* SE Mean */
1248       se_mean = pairs[i].std_dev_diff / sqrt (n) ;
1249       tab_double (trb->t, 4, i+3, TAB_RIGHT, se_mean, NULL);
1250
1251       /* Now work out the confidence interval */
1252       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1253
1254       t = gsl_cdf_tdist_Qinv (q, df);
1255
1256       tab_double (trb->t, 5, i+3, TAB_RIGHT,
1257                 pairs[i].mean_diff - t * se_mean , NULL);
1258
1259       tab_double (trb->t, 6, i+3, TAB_RIGHT,
1260                 pairs[i].mean_diff + t * se_mean , NULL);
1261
1262       t = (pairs[i].mean[0] - pairs[i].mean[1])
1263         / sqrt (
1264                 ( pow2 (pairs[i].s_std_dev[0]) + pow2 (pairs[i].s_std_dev[1]) -
1265                   2 * pairs[i].correlation *
1266                   pairs[i].s_std_dev[0] * pairs[i].s_std_dev[1] )
1267                 / (n - 1)
1268                 );
1269
1270       tab_double (trb->t, 7, i+3, TAB_RIGHT, t, NULL);
1271
1272       /* Degrees of freedom */
1273       tab_double (trb->t, 8, i+3, TAB_RIGHT, df, wfmt);
1274
1275       p = gsl_cdf_tdist_P (t,df);
1276       q = gsl_cdf_tdist_P (t,df);
1277
1278       tab_double (trb->t, 9, i+3, TAB_RIGHT, 2.0* (t>0?q:p), NULL);
1279
1280     }
1281 }
1282
1283 /* Initialize the one sample trbox */
1284 void
1285 trbox_one_sample_init (struct trbox *self, struct cmd_t_test *cmd )
1286 {
1287   const int hsize=7;
1288   const int vsize=cmd->n_variables+3;
1289
1290   self->populate = trbox_one_sample_populate;
1291
1292   trbox_base_init (self, cmd->n_variables,hsize);
1293   tab_title (self->t, _ ("One-Sample Test"));
1294   tab_hline (self->t, TAL_1, 1, hsize - 1, 1);
1295   tab_vline (self->t, TAL_2, 1, 0, vsize - 1);
1296
1297   tab_joint_text (self->t, 1, 0, hsize-1,0, TAB_CENTER | TAT_PRINTF,
1298                  _ ("Test Value = %f"), cmd->n_testval[0]);
1299
1300   tab_box (self->t, -1, -1, -1, TAL_1, 1,1,hsize-1,vsize-1);
1301
1302
1303   tab_joint_text (self->t,5,1,6,1,TAB_CENTER  | TAT_PRINTF,
1304                  _ ("%g%% Confidence Interval of the Difference"),
1305                  cmd->criteria*100.0);
1306
1307   tab_vline (self->t,TAL_GAP,6,1,1);
1308   tab_hline (self->t,TAL_1,5,6,2);
1309   tab_text (self->t, 1, 2, TAB_CENTER | TAT_TITLE, _ ("t"));
1310   tab_text (self->t, 2, 2, TAB_CENTER | TAT_TITLE, _ ("df"));
1311   tab_text (self->t, 3, 2, TAB_CENTER | TAT_TITLE, _ ("Sig. (2-tailed)"));
1312   tab_text (self->t, 4, 2, TAB_CENTER | TAT_TITLE, _ ("Mean Difference"));
1313   tab_text (self->t, 5, 2, TAB_CENTER | TAT_TITLE, _ ("Lower"));
1314   tab_text (self->t, 6, 2, TAB_CENTER | TAT_TITLE, _ ("Upper"));
1315
1316 }
1317
1318
1319 /* Populate the one sample trbox */
1320 static void
1321 trbox_one_sample_populate (struct trbox *trb,
1322                            const struct dictionary *dict,
1323                            struct cmd_t_test *cmd)
1324 {
1325   int i;
1326
1327   const struct variable *wv = dict_get_weight (dict);
1328   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
1329
1330   assert (trb->t);
1331
1332   for (i=0; i < cmd->n_variables; ++i)
1333     {
1334       double t;
1335       double p,q;
1336       double df;
1337       struct group_statistics *gs = &group_proc_get (cmd->v_variables[i])->ugs;
1338
1339
1340       tab_text (trb->t, 0, i+3, TAB_LEFT, var_get_name (cmd->v_variables[i]));
1341
1342       t = (gs->mean - cmd->n_testval[0] ) * sqrt (gs->n) / gs->std_dev ;
1343
1344       tab_double (trb->t, 1, i+3, TAB_RIGHT, t, NULL);
1345
1346       /* degrees of freedom */
1347       df = gs->n - 1;
1348
1349       tab_double (trb->t, 2, i+3, TAB_RIGHT, df, wfmt);
1350
1351       p = gsl_cdf_tdist_P (t, df);
1352       q = gsl_cdf_tdist_Q (t, df);
1353
1354       /* Multiply by 2 to get 2-tailed significance, makeing sure we've got
1355          the correct tail*/
1356       tab_double (trb->t, 3, i+3, TAB_RIGHT, 2.0* (t>0?q:p), NULL);
1357
1358       tab_double (trb->t, 4, i+3, TAB_RIGHT, gs->mean_diff, NULL);
1359
1360
1361       q = (1 - cmd->criteria)/2.0;  /* 2-tailed test */
1362       t = gsl_cdf_tdist_Qinv (q, df);
1363
1364       tab_double (trb->t, 5, i+3, TAB_RIGHT,
1365                  gs->mean_diff - t * gs->se_mean, NULL);
1366
1367       tab_double (trb->t, 6, i+3, TAB_RIGHT,
1368                  gs->mean_diff + t * gs->se_mean, NULL);
1369     }
1370 }
1371
1372 /* Base initializer for the generalized trbox */
1373 void
1374 trbox_base_init (struct trbox *self, size_t data_rows, int cols)
1375 {
1376   const size_t rows = 3 + data_rows;
1377
1378   self->finalize = trbox_base_finalize;
1379   self->t = tab_create (cols, rows, 0);
1380   tab_headers (self->t,0,0,3,0);
1381   tab_box (self->t, TAL_2, TAL_2, TAL_0, TAL_0, 0, 0, cols -1, rows -1);
1382   tab_hline (self->t, TAL_2,0,cols-1,3);
1383   tab_dim (self->t, tab_natural_dimensions);
1384 }
1385
1386
1387 /* Base finalizer for the trbox */
1388 void
1389 trbox_base_finalize (struct trbox *trb)
1390 {
1391   tab_submit (trb->t);
1392 }
1393
1394
1395 /* Create , populate and submit the Paired Samples Correlation box */
1396 static void
1397 pscbox (const struct dictionary *dict)
1398 {
1399   const struct variable *wv = dict_get_weight (dict);
1400   const struct fmt_spec *wfmt = wv ? var_get_print_format (wv) : &F_8_0;
1401
1402   const int rows = 1 + n_pairs;
1403   const int cols = 5;
1404   int i;
1405
1406   struct tab_table *table;
1407
1408   table = tab_create (cols,rows,0);
1409
1410   tab_columns (table, SOM_COL_DOWN, 1);
1411   tab_headers (table,0,0,1,0);
1412   tab_box (table, TAL_2, TAL_2, TAL_0, TAL_1, 0, 0, cols -1, rows -1 );
1413   tab_hline (table, TAL_2, 0, cols - 1, 1);
1414   tab_vline (table, TAL_2, 2, 0, rows - 1);
1415   tab_dim (table, tab_natural_dimensions);
1416   tab_title (table, _ ("Paired Samples Correlations"));
1417
1418   /* column headings */
1419   tab_text (table, 2,0, TAB_CENTER | TAT_TITLE, _ ("N"));
1420   tab_text (table, 3,0, TAB_CENTER | TAT_TITLE, _ ("Correlation"));
1421   tab_text (table, 4,0, TAB_CENTER | TAT_TITLE, _ ("Sig."));
1422
1423   for (i=0; i < n_pairs; ++i)
1424     {
1425       double p,q;
1426
1427       double df = pairs[i].n -2;
1428
1429       double correlation_t =
1430         pairs[i].correlation * sqrt (df) /
1431         sqrt (1 - pow2 (pairs[i].correlation));
1432
1433
1434       /* row headings */
1435       tab_text (table, 0,i+1, TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1436                _ ("Pair %d"), i);
1437
1438       tab_text (table, 1,i+1, TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1439                _ ("%s & %s"),
1440                var_get_name (pairs[i].v[0]),
1441                var_get_name (pairs[i].v[1]));
1442
1443
1444       /* row data */
1445       tab_double (table, 2, i+1, TAB_RIGHT, pairs[i].n, wfmt);
1446       tab_double (table, 3, i+1, TAB_RIGHT, pairs[i].correlation, NULL);
1447
1448       p = gsl_cdf_tdist_P (correlation_t, df);
1449       q = gsl_cdf_tdist_Q (correlation_t, df);
1450
1451       tab_double (table, 4, i+1, TAB_RIGHT, 2.0* (correlation_t>0?q:p), NULL);
1452     }
1453
1454   tab_submit (table);
1455 }
1456
1457
1458
1459
1460 /* Calculation Implementation */
1461
1462 /* Per case calculations common to all variants of the T test */
1463 static int
1464 common_calc (const struct dictionary *dict,
1465              const struct ccase *c,
1466              void *_cmd,
1467              enum mv_class exclude)
1468 {
1469   int i;
1470   struct cmd_t_test *cmd = (struct cmd_t_test *)_cmd;
1471
1472   double weight = dict_get_case_weight (dict, c, NULL);
1473
1474
1475   /* Listwise has to be implicit if the independent variable is missing ?? */
1476   if ( cmd->sbc_groups )
1477     {
1478       if (var_is_value_missing (indep_var, case_data (c, indep_var), exclude))
1479         return 0;
1480     }
1481
1482   for (i = 0; i < cmd->n_variables ; ++i)
1483     {
1484       const struct variable *v = cmd->v_variables[i];
1485       const union value *val = case_data (c, v);
1486
1487       if (!var_is_value_missing (v, val, exclude))
1488         {
1489           struct group_statistics *gs;
1490           gs = &group_proc_get (v)->ugs;
1491
1492           gs->n += weight;
1493           gs->sum += weight * val->f;
1494           gs->ssq += weight * pow2 (val->f);
1495         }
1496     }
1497   return 0;
1498 }
1499
1500 /* Pre calculations common to all variants of the T test */
1501 static void
1502 common_precalc ( struct cmd_t_test *cmd )
1503 {
1504   int i=0;
1505
1506   for (i=0; i< cmd->n_variables ; ++i)
1507     {
1508       struct group_statistics *gs;
1509       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1510
1511       gs->sum=0;
1512       gs->n=0;
1513       gs->ssq=0;
1514       gs->sum_diff=0;
1515     }
1516 }
1517
1518 /* Post calculations common to all variants of the T test */
1519 void
1520 common_postcalc (struct cmd_t_test *cmd)
1521 {
1522   int i=0;
1523
1524   for (i=0; i< cmd->n_variables ; ++i)
1525     {
1526       struct group_statistics *gs;
1527       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1528
1529       gs->mean=gs->sum / gs->n;
1530       gs->s_std_dev= sqrt (
1531                          ( (gs->ssq / gs->n ) - pow2 (gs->mean))
1532                          ) ;
1533
1534       gs->std_dev= sqrt (
1535                          gs->n/ (gs->n-1) *
1536                          ( (gs->ssq / gs->n ) - pow2 (gs->mean))
1537                          ) ;
1538
1539       gs->se_mean = gs->std_dev / sqrt (gs->n);
1540       gs->mean_diff= gs->sum_diff / gs->n;
1541     }
1542 }
1543
1544 /* Per case calculations for one sample t test  */
1545 static int
1546 one_sample_calc (const struct dictionary *dict,
1547                  const struct ccase *c, void *cmd_,
1548                  enum mv_class exclude)
1549 {
1550   int i;
1551
1552   struct cmd_t_test *cmd = (struct cmd_t_test *)cmd_;
1553
1554   double weight = dict_get_case_weight (dict, c, NULL);
1555
1556
1557   for (i=0; i< cmd->n_variables ; ++i)
1558     {
1559       struct group_statistics *gs;
1560       const struct variable *v = cmd->v_variables[i];
1561       const union value *val = case_data (c, v);
1562
1563       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1564
1565       if (!var_is_value_missing (v, val, exclude))
1566         gs->sum_diff += weight * (val->f - cmd->n_testval[0]);
1567     }
1568
1569   return 0;
1570 }
1571
1572 /* Pre calculations for one sample t test */
1573 static void
1574 one_sample_precalc ( struct cmd_t_test *cmd )
1575 {
1576   int i=0;
1577
1578   for (i=0; i< cmd->n_variables ; ++i)
1579     {
1580       struct group_statistics *gs;
1581       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1582
1583       gs->sum_diff=0;
1584     }
1585 }
1586
1587 /* Post calculations for one sample t test */
1588 static void
1589 one_sample_postcalc (struct cmd_t_test *cmd)
1590 {
1591   int i=0;
1592
1593   for (i=0; i< cmd->n_variables ; ++i)
1594     {
1595       struct group_statistics *gs;
1596       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1597
1598       gs->mean_diff = gs->sum_diff / gs->n ;
1599     }
1600 }
1601
1602
1603
1604 static void
1605 paired_precalc (struct cmd_t_test *cmd UNUSED)
1606 {
1607   int i;
1608
1609   for (i=0; i < n_pairs ; ++i )
1610     {
1611       pairs[i].n = 0;
1612       pairs[i].sum[0] = 0;      pairs[i].sum[1] = 0;
1613       pairs[i].ssq[0] = 0;      pairs[i].ssq[1] = 0;
1614       pairs[i].sum_of_prod = 0;
1615       pairs[i].correlation = 0;
1616       pairs[i].sum_of_diffs = 0;
1617       pairs[i].ssq_diffs = 0;
1618     }
1619
1620 }
1621
1622
1623 static int
1624 paired_calc (const struct dictionary *dict, const struct ccase *c,
1625              struct cmd_t_test *cmd UNUSED, enum mv_class exclude)
1626 {
1627   int i;
1628
1629   double weight = dict_get_case_weight (dict, c, NULL);
1630
1631   for (i=0; i < n_pairs ; ++i )
1632     {
1633       const struct variable *v0 = pairs[i].v[0];
1634       const struct variable *v1 = pairs[i].v[1];
1635
1636       const union value *val0 = case_data (c, v0);
1637       const union value *val1 = case_data (c, v1);
1638
1639       if (!var_is_value_missing (v0, val0, exclude) &&
1640           !var_is_value_missing (v1, val1, exclude))
1641         {
1642           pairs[i].n += weight;
1643           pairs[i].sum[0] += weight * val0->f;
1644           pairs[i].sum[1] += weight * val1->f;
1645
1646           pairs[i].ssq[0] += weight * pow2 (val0->f);
1647           pairs[i].ssq[1] += weight * pow2 (val1->f);
1648
1649           pairs[i].sum_of_prod += weight * val0->f * val1->f ;
1650
1651           pairs[i].sum_of_diffs += weight * ( val0->f - val1->f ) ;
1652           pairs[i].ssq_diffs += weight * pow2 (val0->f - val1->f);
1653         }
1654     }
1655
1656   return 0;
1657 }
1658
1659 static void
1660 paired_postcalc (struct cmd_t_test *cmd UNUSED)
1661 {
1662   int i;
1663
1664   for (i=0; i < n_pairs ; ++i )
1665     {
1666       int j;
1667       const double n = pairs[i].n;
1668
1669       for (j=0; j < 2 ; ++j)
1670         {
1671           pairs[i].mean[j] = pairs[i].sum[j] / n ;
1672           pairs[i].s_std_dev[j] = sqrt ((pairs[i].ssq[j] / n -
1673                                               pow2 (pairs[i].mean[j]))
1674                                      );
1675
1676           pairs[i].std_dev[j] = sqrt (n/ (n-1)* (pairs[i].ssq[j] / n -
1677                                               pow2 (pairs[i].mean[j]))
1678                                      );
1679         }
1680
1681       pairs[i].correlation = pairs[i].sum_of_prod / pairs[i].n -
1682         pairs[i].mean[0] * pairs[i].mean[1] ;
1683       /* correlation now actually contains the covariance */
1684
1685       pairs[i].correlation /= pairs[i].std_dev[0] * pairs[i].std_dev[1];
1686       pairs[i].correlation *= pairs[i].n / ( pairs[i].n - 1 );
1687
1688       pairs[i].mean_diff = pairs[i].sum_of_diffs / n ;
1689
1690       pairs[i].std_dev_diff = sqrt (  n / (n - 1) * (
1691                                     ( pairs[i].ssq_diffs / n )
1692                                     -
1693                                     pow2 (pairs[i].mean_diff )
1694                                     ) );
1695     }
1696 }
1697
1698 static void
1699 group_precalc (struct cmd_t_test *cmd )
1700 {
1701   int i;
1702   int j;
1703
1704   for (i=0; i< cmd->n_variables ; ++i)
1705     {
1706       struct group_proc *ttpr = group_proc_get (cmd->v_variables[i]);
1707
1708       /* There's always 2 groups for a T - TEST */
1709       ttpr->n_groups = 2;
1710
1711       gp.indep_var = indep_var;
1712
1713       ttpr->group_hash = hsh_create (2,
1714                                     (hsh_compare_func *) compare_group_binary,
1715                                     (hsh_hash_func *) hash_group_binary,
1716                                     (hsh_free_func *) free_group,
1717                                     (void *) &gp );
1718
1719       for (j=0 ; j < 2 ; ++j)
1720         {
1721           struct group_statistics *gs = xmalloc (sizeof *gs);
1722
1723           gs->sum = 0;
1724           gs->n = 0;
1725           gs->ssq = 0;
1726
1727           if ( gp.criterion == CMP_EQ )
1728             {
1729               gs->id = gp.v.g_value[j];
1730             }
1731           else
1732             {
1733               if ( j == 0 )
1734                 gs->id.f = gp.v.critical_value - 1.0 ;
1735               else
1736                 gs->id.f = gp.v.critical_value + 1.0 ;
1737             }
1738
1739           hsh_insert ( ttpr->group_hash, (void *) gs );
1740         }
1741     }
1742
1743 }
1744
1745 static int
1746 group_calc (const struct dictionary *dict,
1747             const struct ccase *c, struct cmd_t_test *cmd,
1748             enum mv_class exclude)
1749 {
1750   int i;
1751
1752   const double weight = dict_get_case_weight (dict, c, NULL);
1753
1754   const union value *gv;
1755
1756   if (var_is_value_missing (indep_var, case_data (c, indep_var), exclude))
1757     return 0;
1758
1759   gv = case_data (c, indep_var);
1760
1761   for (i=0; i< cmd->n_variables ; ++i)
1762     {
1763       const struct variable *var = cmd->v_variables[i];
1764       const union value *val = case_data (c, var);
1765       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1766       struct group_statistics *gs;
1767
1768       gs = hsh_find (grp_hash, (void *) gv);
1769
1770       /* If the independent variable doesn't match either of the values
1771          for this case then move on to the next case */
1772       if ( ! gs )
1773         return 0;
1774
1775       if (!var_is_value_missing (var, val, exclude))
1776         {
1777           gs->n += weight;
1778           gs->sum += weight * val->f;
1779           gs->ssq += weight * pow2 (val->f);
1780         }
1781     }
1782
1783   return 0;
1784 }
1785
1786
1787 static void
1788 group_postcalc ( struct cmd_t_test *cmd )
1789 {
1790   int i;
1791
1792   for (i = 0; i < cmd->n_variables ; ++i)
1793     {
1794       const struct variable *var = cmd->v_variables[i];
1795       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1796       struct hsh_iterator g;
1797       struct group_statistics *gs;
1798       int count=0;
1799
1800       for (gs =  hsh_first (grp_hash,&g);
1801            gs != 0;
1802            gs = hsh_next (grp_hash,&g))
1803         {
1804           gs->mean = gs->sum / gs->n;
1805
1806           gs->s_std_dev= sqrt (
1807                               ( (gs->ssq / gs->n ) - pow2 (gs->mean))
1808                               ) ;
1809
1810           gs->std_dev= sqrt (
1811                             gs->n/ (gs->n-1) *
1812                             ( (gs->ssq / gs->n ) - pow2 (gs->mean))
1813                             ) ;
1814
1815           gs->se_mean = gs->std_dev / sqrt (gs->n);
1816           count ++;
1817         }
1818       assert (count == 2);
1819     }
1820 }
1821
1822
1823
1824 static void
1825 calculate (struct cmd_t_test *cmd,
1826           struct casereader *input, const struct dataset *ds)
1827 {
1828   const struct dictionary *dict = dataset_dict (ds);
1829   struct ssbox stat_summary_box;
1830   struct trbox test_results_box;
1831
1832   struct casereader *pass1, *pass2, *pass3;
1833   struct taint *taint;
1834   struct ccase *c;
1835
1836   enum mv_class exclude = cmd->miss != TTS_INCLUDE ? MV_ANY : MV_SYSTEM;
1837
1838   c = casereader_peek (input, 0);
1839   if (c == NULL)
1840     {
1841       casereader_destroy (input);
1842       return;
1843     }
1844   output_split_file_values (ds, c);
1845   case_unref (c);
1846
1847   if ( cmd->miss == TTS_LISTWISE )
1848     input = casereader_create_filter_missing (input,
1849                                               cmd->v_variables,
1850                                               cmd->n_variables,
1851                                               exclude, NULL, NULL);
1852
1853   input = casereader_create_filter_weight (input, dict, NULL, NULL);
1854
1855   taint = taint_clone (casereader_get_taint (input));
1856   casereader_split (input, &pass1, &pass2);
1857
1858   common_precalc (cmd);
1859   for (; (c = casereader_read (pass1)) != NULL; case_unref (c))
1860     common_calc (dict, c, cmd, exclude);
1861   casereader_destroy (pass1);
1862   common_postcalc (cmd);
1863
1864   switch (mode)
1865     {
1866     case T_1_SAMPLE:
1867       one_sample_precalc (cmd);
1868       for (; (c = casereader_read (pass2)) != NULL; case_unref (c))
1869         one_sample_calc (dict, c, cmd, exclude);
1870       one_sample_postcalc (cmd);
1871       break;
1872     case T_PAIRED:
1873       paired_precalc (cmd);
1874       for (; (c = casereader_read (pass2)) != NULL; case_unref (c))
1875         paired_calc (dict, c, cmd, exclude);
1876       paired_postcalc (cmd);
1877       break;
1878     case T_IND_SAMPLES:
1879       pass3 = casereader_clone (pass2);
1880
1881       group_precalc (cmd);
1882       for (; (c = casereader_read (pass2)) != NULL; case_unref (c))
1883         group_calc (dict, c, cmd, exclude);
1884       group_postcalc (cmd);
1885
1886       levene (dict, pass3, indep_var, cmd->n_variables, cmd->v_variables,
1887               exclude);
1888       break;
1889     }
1890   casereader_destroy (pass2);
1891
1892   if (!taint_has_tainted_successor (taint))
1893     {
1894       ssbox_create (&stat_summary_box,cmd,mode);
1895       ssbox_populate (&stat_summary_box, dict, cmd);
1896       ssbox_finalize (&stat_summary_box);
1897
1898       if ( mode == T_PAIRED )
1899         pscbox (dict);
1900
1901       trbox_create (&test_results_box, cmd, mode);
1902       trbox_populate (&test_results_box, dict, cmd);
1903       trbox_finalize (&test_results_box);
1904     }
1905
1906   taint_destroy (taint);
1907 }
1908
1909 short which_group (const struct group_statistics *g,
1910                   const struct group_properties *p);
1911
1912 /* Return -1 if the id of a is less than b; +1 if greater than and
1913    0 if equal */
1914 static int
1915 compare_group_binary (const struct group_statistics *a,
1916                      const struct group_statistics *b,
1917                      const struct group_properties *p)
1918 {
1919   short flag_a;
1920   short flag_b;
1921
1922   if ( p->criterion == CMP_LE )
1923     {
1924       flag_a = ( a->id.f < p->v.critical_value ) ;
1925       flag_b = ( b->id.f < p->v.critical_value ) ;
1926     }
1927   else
1928     {
1929       flag_a = which_group (a, p);
1930       flag_b = which_group (b, p);
1931     }
1932
1933   if (flag_a < flag_b )
1934     return -1;
1935
1936   return (flag_a > flag_b);
1937 }
1938
1939 /* This is a degenerate case of a hash, since it can only return three possible
1940    values.  It's really a comparison, being used as a hash function */
1941
1942 static unsigned
1943 hash_group_binary (const struct group_statistics *g,
1944                   const struct group_properties *p)
1945 {
1946   short flag = -1;
1947
1948   if ( p->criterion == CMP_LE )
1949     {
1950       flag = ( g->id.f < p->v.critical_value ) ;
1951     }
1952   else if ( p->criterion == CMP_EQ)
1953     {
1954       flag = which_group (g,p);
1955     }
1956   else
1957     NOT_REACHED ();
1958
1959   return flag;
1960 }
1961
1962 /* return 0 if G belongs to group 0,
1963           1 if it belongs to group 1,
1964           2 if it belongs to neither group */
1965 short
1966 which_group (const struct group_statistics *g,
1967             const struct group_properties *p)
1968 {
1969   if ( 0 == compare_values_short (&g->id, &p->v.g_value[0], p->indep_var))
1970     return 0;
1971
1972   if ( 0 == compare_values_short (&g->id, &p->v.g_value[1], p->indep_var))
1973     return 1;
1974
1975   return 2;
1976 }
1977
1978 /*
1979   Local Variables:
1980   mode: c
1981   End:
1982 */