pintos-os.org Git - pspp/blob - 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 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_Q (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       /* corr2 will mathematically always be in the range [0, 1.0].  Inaccurate
1436          calculations sometimes cause it to be slightly greater than 1.0, so
1437          force it into the correct range to avoid NaN from sqrt(). */
1438       double corr2 = MIN (1.0, pow2 (pairs[i].correlation));
1439       double correlation_t =
1440         pairs[i].correlation * sqrt (df) /
1441         sqrt (1 - corr2);
1442
1443
1444       /* row headings */
1445       tab_text (table, 0,i+1, TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1446                _ ("Pair %d"), i);
1447
1448       tab_text (table, 1,i+1, TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1449                _ ("%s & %s"),
1450                var_get_name (pairs[i].v[0]),
1451                var_get_name (pairs[i].v[1]));
1452
1453
1454       /* row data */
1455       tab_double (table, 2, i+1, TAB_RIGHT, pairs[i].n, wfmt);
1456       tab_double (table, 3, i+1, TAB_RIGHT, pairs[i].correlation, NULL);
1457
1458       p = gsl_cdf_tdist_P (correlation_t, df);
1459       q = gsl_cdf_tdist_Q (correlation_t, df);
1460
1461       tab_double (table, 4, i+1, TAB_RIGHT, 2.0* (correlation_t>0?q:p), NULL);
1462     }
1463
1464   tab_submit (table);
1465 }
1466
1467
1468
1469
1470 /* Calculation Implementation */
1471
1472 /* Per case calculations common to all variants of the T test */
1473 static int
1474 common_calc (const struct dictionary *dict,
1475              const struct ccase *c,
1476              void *_cmd,
1477              enum mv_class exclude)
1478 {
1479   int i;
1480   struct cmd_t_test *cmd = (struct cmd_t_test *)_cmd;
1481
1482   double weight = dict_get_case_weight (dict, c, NULL);
1483
1484
1485   /* Listwise has to be implicit if the independent variable is missing ?? */
1486   if ( cmd->sbc_groups )
1487     {
1488       if (var_is_value_missing (indep_var, case_data (c, indep_var), exclude))
1489         return 0;
1490     }
1491
1492   for (i = 0; i < cmd->n_variables ; ++i)
1493     {
1494       const struct variable *v = cmd->v_variables[i];
1495       const union value *val = case_data (c, v);
1496
1497       if (!var_is_value_missing (v, val, exclude))
1498         {
1499           struct group_statistics *gs;
1500           gs = &group_proc_get (v)->ugs;
1501
1502           gs->n += weight;
1503           gs->sum += weight * val->f;
1504           gs->ssq += weight * val->f * val->f;
1505         }
1506     }
1507   return 0;
1508 }
1509
1510 /* Pre calculations common to all variants of the T test */
1511 static void
1512 common_precalc ( struct cmd_t_test *cmd )
1513 {
1514   int i=0;
1515
1516   for (i=0; i< cmd->n_variables ; ++i)
1517     {
1518       struct group_statistics *gs;
1519       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1520
1521       gs->sum=0;
1522       gs->n=0;
1523       gs->ssq=0;
1524       gs->sum_diff=0;
1525     }
1526 }
1527
1528 /* Post calculations common to all variants of the T test */
1529 void
1530 common_postcalc (struct cmd_t_test *cmd)
1531 {
1532   int i=0;
1533
1534   for (i=0; i< cmd->n_variables ; ++i)
1535     {
1536       struct group_statistics *gs;
1537       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1538
1539       gs->mean=gs->sum / gs->n;
1540       gs->s_std_dev= sqrt (
1541                          ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1542                          ) ;
1543
1544       gs->std_dev= sqrt (
1545                          gs->n/ (gs->n-1) *
1546                          ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1547                          ) ;
1548
1549       gs->se_mean = gs->std_dev / sqrt (gs->n);
1550       gs->mean_diff= gs->sum_diff / gs->n;
1551     }
1552 }
1553
1554 /* Per case calculations for one sample t test  */
1555 static int
1556 one_sample_calc (const struct dictionary *dict,
1557                  const struct ccase *c, void *cmd_,
1558                  enum mv_class exclude)
1559 {
1560   int i;
1561
1562   struct cmd_t_test *cmd = (struct cmd_t_test *)cmd_;
1563
1564   double weight = dict_get_case_weight (dict, c, NULL);
1565
1566
1567   for (i=0; i< cmd->n_variables ; ++i)
1568     {
1569       struct group_statistics *gs;
1570       const struct variable *v = cmd->v_variables[i];
1571       const union value *val = case_data (c, v);
1572
1573       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1574
1575       if (!var_is_value_missing (v, val, exclude))
1576         gs->sum_diff += weight * (val->f - cmd->n_testval[0]);
1577     }
1578
1579   return 0;
1580 }
1581
1582 /* Pre calculations for one sample t test */
1583 static void
1584 one_sample_precalc ( struct cmd_t_test *cmd )
1585 {
1586   int i=0;
1587
1588   for (i=0; i< cmd->n_variables ; ++i)
1589     {
1590       struct group_statistics *gs;
1591       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1592
1593       gs->sum_diff=0;
1594     }
1595 }
1596
1597 /* Post calculations for one sample t test */
1598 static void
1599 one_sample_postcalc (struct cmd_t_test *cmd)
1600 {
1601   int i=0;
1602
1603   for (i=0; i< cmd->n_variables ; ++i)
1604     {
1605       struct group_statistics *gs;
1606       gs= &group_proc_get (cmd->v_variables[i])->ugs;
1607
1608       gs->mean_diff = gs->sum_diff / gs->n ;
1609     }
1610 }
1611
1612
1613
1614 static void
1615 paired_precalc (struct cmd_t_test *cmd UNUSED)
1616 {
1617   int i;
1618
1619   for (i=0; i < n_pairs ; ++i )
1620     {
1621       pairs[i].n = 0;
1622       pairs[i].sum[0] = 0;      pairs[i].sum[1] = 0;
1623       pairs[i].ssq[0] = 0;      pairs[i].ssq[1] = 0;
1624       pairs[i].sum_of_prod = 0;
1625       pairs[i].correlation = 0;
1626       pairs[i].sum_of_diffs = 0;
1627       pairs[i].ssq_diffs = 0;
1628     }
1629
1630 }
1631
1632
1633 static int
1634 paired_calc (const struct dictionary *dict, const struct ccase *c,
1635              struct cmd_t_test *cmd UNUSED, enum mv_class exclude)
1636 {
1637   int i;
1638
1639   double weight = dict_get_case_weight (dict, c, NULL);
1640
1641   for (i=0; i < n_pairs ; ++i )
1642     {
1643       const struct variable *v0 = pairs[i].v[0];
1644       const struct variable *v1 = pairs[i].v[1];
1645
1646       const union value *val0 = case_data (c, v0);
1647       const union value *val1 = case_data (c, v1);
1648
1649       if (!var_is_value_missing (v0, val0, exclude) &&
1650           !var_is_value_missing (v1, val1, exclude))
1651         {
1652           pairs[i].n += weight;
1653           pairs[i].sum[0] += weight * val0->f;
1654           pairs[i].sum[1] += weight * val1->f;
1655
1656           pairs[i].ssq[0] += weight * pow2 (val0->f);
1657           pairs[i].ssq[1] += weight * pow2 (val1->f);
1658
1659           pairs[i].sum_of_prod += weight * val0->f * val1->f ;
1660
1661           pairs[i].sum_of_diffs += weight * ( val0->f - val1->f ) ;
1662           pairs[i].ssq_diffs += weight * pow2 (val0->f - val1->f);
1663         }
1664     }
1665
1666   return 0;
1667 }
1668
1669 static void
1670 paired_postcalc (struct cmd_t_test *cmd UNUSED)
1671 {
1672   int i;
1673
1674   for (i=0; i < n_pairs ; ++i )
1675     {
1676       int j;
1677       const double n = pairs[i].n;
1678
1679       for (j=0; j < 2 ; ++j)
1680         {
1681           pairs[i].mean[j] = pairs[i].sum[j] / n ;
1682           pairs[i].s_std_dev[j] = sqrt ((pairs[i].ssq[j] / n -
1683                                               pow2 (pairs[i].mean[j]))
1684                                      );
1685
1686           pairs[i].std_dev[j] = sqrt (n/ (n-1)* (pairs[i].ssq[j] / n -
1687                                               pow2 (pairs[i].mean[j]))
1688                                      );
1689         }
1690
1691       pairs[i].correlation = pairs[i].sum_of_prod / pairs[i].n -
1692         pairs[i].mean[0] * pairs[i].mean[1] ;
1693       /* correlation now actually contains the covariance */
1694
1695       pairs[i].correlation /= pairs[i].std_dev[0] * pairs[i].std_dev[1];
1696       pairs[i].correlation *= pairs[i].n / ( pairs[i].n - 1 );
1697
1698       pairs[i].mean_diff = pairs[i].sum_of_diffs / n ;
1699
1700       pairs[i].std_dev_diff = sqrt (  n / (n - 1) * (
1701                                     ( pairs[i].ssq_diffs / n )
1702                                     -
1703                                     pow2 (pairs[i].mean_diff )
1704                                     ) );
1705     }
1706 }
1707
1708 static void
1709 group_precalc (struct cmd_t_test *cmd )
1710 {
1711   int i;
1712   int j;
1713
1714   for (i=0; i< cmd->n_variables ; ++i)
1715     {
1716       struct group_proc *ttpr = group_proc_get (cmd->v_variables[i]);
1717
1718       /* There's always 2 groups for a T - TEST */
1719       ttpr->n_groups = 2;
1720
1721       gp.indep_width = var_get_width (indep_var);
1722
1723       ttpr->group_hash = hsh_create (2,
1724                                     (hsh_compare_func *) compare_group_binary,
1725                                     (hsh_hash_func *) hash_group_binary,
1726                                     (hsh_free_func *) free_group,
1727                                     (void *) &gp );
1728
1729       for (j=0 ; j < 2 ; ++j)
1730         {
1731           struct group_statistics *gs = xmalloc (sizeof *gs);
1732
1733           gs->sum = 0;
1734           gs->n = 0;
1735           gs->ssq = 0;
1736
1737           if ( gp.criterion == CMP_EQ )
1738             {
1739               gs->id = gp.v.g_value[j];
1740             }
1741           else
1742             {
1743               if ( j == 0 )
1744                 gs->id.f = gp.v.critical_value - 1.0 ;
1745               else
1746                 gs->id.f = gp.v.critical_value + 1.0 ;
1747             }
1748
1749           hsh_insert ( ttpr->group_hash, (void *) gs );
1750         }
1751     }
1752
1753 }
1754
1755 static int
1756 group_calc (const struct dictionary *dict,
1757             const struct ccase *c, struct cmd_t_test *cmd,
1758             enum mv_class exclude)
1759 {
1760   int i;
1761
1762   const double weight = dict_get_case_weight (dict, c, NULL);
1763
1764   const union value *gv;
1765
1766   if (var_is_value_missing (indep_var, case_data (c, indep_var), exclude))
1767     return 0;
1768
1769   gv = case_data (c, indep_var);
1770
1771   for (i=0; i< cmd->n_variables ; ++i)
1772     {
1773       const struct variable *var = cmd->v_variables[i];
1774       const union value *val = case_data (c, var);
1775       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1776       struct group_statistics *gs;
1777
1778       gs = hsh_find (grp_hash, (void *) gv);
1779
1780       /* If the independent variable doesn't match either of the values
1781          for this case then move on to the next case */
1782       if ( ! gs )
1783         return 0;
1784
1785       if (!var_is_value_missing (var, val, exclude))
1786         {
1787           gs->n += weight;
1788           gs->sum += weight * val->f;
1789           gs->ssq += weight * pow2 (val->f);
1790         }
1791     }
1792
1793   return 0;
1794 }
1795
1796
1797 static void
1798 group_postcalc ( struct cmd_t_test *cmd )
1799 {
1800   int i;
1801
1802   for (i = 0; i < cmd->n_variables ; ++i)
1803     {
1804       const struct variable *var = cmd->v_variables[i];
1805       struct hsh_table *grp_hash = group_proc_get (var)->group_hash;
1806       struct hsh_iterator g;
1807       struct group_statistics *gs;
1808       int count=0;
1809
1810       for (gs =  hsh_first (grp_hash,&g);
1811            gs != 0;
1812            gs = hsh_next (grp_hash,&g))
1813         {
1814           gs->mean = gs->sum / gs->n;
1815
1816           gs->s_std_dev= sqrt (
1817                               ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1818                               ) ;
1819
1820           gs->std_dev= sqrt (
1821                             gs->n/ (gs->n-1) *
1822                             ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1823                             ) ;
1824
1825           gs->se_mean = gs->std_dev / sqrt (gs->n);
1826           count ++;
1827         }
1828       assert (count == 2);
1829     }
1830 }
1831
1832
1833
1834 static void
1835 calculate (struct cmd_t_test *cmd,
1836           struct casereader *input, const struct dataset *ds)
1837 {
1838   const struct dictionary *dict = dataset_dict (ds);
1839   struct ssbox stat_summary_box;
1840   struct trbox test_results_box;
1841
1842   struct casereader *pass1, *pass2, *pass3;
1843   struct taint *taint;
1844   struct ccase c;
1845
1846   enum mv_class exclude = cmd->miss != TTS_INCLUDE ? MV_ANY : MV_SYSTEM;
1847
1848   if (!casereader_peek (input, 0, &c))
1849     {
1850       casereader_destroy (input);
1851       return;
1852     }
1853   output_split_file_values (ds, &c);
1854   case_destroy (&c);
1855
1856   if ( cmd->miss == TTS_LISTWISE )
1857     input = casereader_create_filter_missing (input,
1858                                               cmd->v_variables,
1859                                               cmd->n_variables,
1860                                               exclude, NULL);
1861
1862   input = casereader_create_filter_weight (input, dict, NULL, NULL);
1863
1864   taint = taint_clone (casereader_get_taint (input));
1865   casereader_split (input, &pass1, &pass2);
1866
1867   common_precalc (cmd);
1868   for (; casereader_read (pass1, &c); case_destroy (&c))
1869     common_calc (dict, &c, cmd, exclude);
1870   casereader_destroy (pass1);
1871   common_postcalc (cmd);
1872
1873   switch (mode)
1874     {
1875     case T_1_SAMPLE:
1876       one_sample_precalc (cmd);
1877       for (; casereader_read (pass2, &c); case_destroy (&c))
1878         one_sample_calc (dict, &c, cmd, exclude);
1879       one_sample_postcalc (cmd);
1880       break;
1881     case T_PAIRED:
1882       paired_precalc (cmd);
1883       for (; casereader_read (pass2, &c); case_destroy (&c))
1884         paired_calc (dict, &c, cmd, exclude);
1885       paired_postcalc (cmd);
1886       break;
1887     case T_IND_SAMPLES:
1888       pass3 = casereader_clone (pass2);
1889
1890       group_precalc (cmd);
1891       for (; casereader_read (pass2, &c); case_destroy (&c))
1892         group_calc (dict, &c, cmd, exclude);
1893       group_postcalc (cmd);
1894
1895       levene (dict, pass3, indep_var, cmd->n_variables, cmd->v_variables,
1896               exclude);
1897       break;
1898     }
1899   casereader_destroy (pass2);
1900
1901   if (!taint_has_tainted_successor (taint))
1902     {
1903       ssbox_create (&stat_summary_box,cmd,mode);
1904       ssbox_populate (&stat_summary_box, dict, cmd);
1905       ssbox_finalize (&stat_summary_box);
1906
1907       if ( mode == T_PAIRED )
1908         pscbox (dict);
1909
1910       trbox_create (&test_results_box, cmd, mode);
1911       trbox_populate (&test_results_box, dict, cmd);
1912       trbox_finalize (&test_results_box);
1913     }
1914
1915   taint_destroy (taint);
1916 }
1917
1918 short which_group (const struct group_statistics *g,
1919                   const struct group_properties *p);
1920
1921 /* Return -1 if the id of a is less than b; +1 if greater than and
1922    0 if equal */
1923 static int
1924 compare_group_binary (const struct group_statistics *a,
1925                      const struct group_statistics *b,
1926                      const struct group_properties *p)
1927 {
1928   short flag_a;
1929   short flag_b;
1930
1931   if ( p->criterion == CMP_LE )
1932     {
1933       /* less-than comparision is not meaningfull for
1934          alpha variables, so we shouldn't ever arrive here */
1935       assert (p->indep_width == 0 ) ;
1936
1937       flag_a = ( a->id.f < p->v.critical_value ) ;
1938       flag_b = ( b->id.f < p->v.critical_value ) ;
1939     }
1940   else
1941     {
1942       flag_a = which_group (a, p);
1943       flag_b = which_group (b, p);
1944     }
1945
1946   if (flag_a < flag_b )
1947     return -1;
1948
1949   return (flag_a > flag_b);
1950 }
1951
1952 /* This is a degenerate case of a hash, since it can only return three possible
1953    values.  It's really a comparison, being used as a hash function */
1954
1955 static unsigned
1956 hash_group_binary (const struct group_statistics *g,
1957                   const struct group_properties *p)
1958 {
1959   short flag = -1;
1960
1961   if ( p->criterion == CMP_LE )
1962     {
1963       /* Not meaningfull to do a less than compare for alpha values ? */
1964       assert (p->indep_width == 0 ) ;
1965       flag = ( g->id.f < p->v.critical_value ) ;
1966     }
1967   else if ( p->criterion == CMP_EQ)
1968     {
1969       flag = which_group (g,p);
1970     }
1971   else
1972     NOT_REACHED ();
1973
1974   return flag;
1975 }
1976
1977 /* return 0 if G belongs to group 0,
1978           1 if it belongs to group 1,
1979           2 if it belongs to neither group */
1980 short
1981 which_group (const struct group_statistics *g,
1982             const struct group_properties *p)
1983 {
1984   if ( 0 == compare_values (&g->id, &p->v.g_value[0], p->indep_width))
1985     return 0;
1986
1987   if ( 0 == compare_values (&g->id, &p->v.g_value[1], p->indep_width))
1988     return 1;
1989
1990   return 2;
1991 }
1992
1993 /*
1994   Local Variables:
1995   mode: c
1996   End:
1997 */