src/language/stats/t-test.q

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