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