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