f33a990edf06797ec401fa94116eeaf696f37de4
[pspp-builds.git] / src / language / stats / oneway.q
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 1997-9, 2000, 2007 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 <stdio.h>
22 #include <stdlib.h>
23
24 #include <data/case.h>
25 #include <data/casegrouper.h>
26 #include <data/casereader.h>
27 #include <data/dictionary.h>
28 #include <data/procedure.h>
29 #include <data/value-labels.h>
30 #include <data/variable.h>
31 #include <language/command.h>
32 #include <language/dictionary/split-file.h>
33 #include <language/lexer/lexer.h>
34 #include <libpspp/alloc.h>
35 #include <libpspp/compiler.h>
36 #include <libpspp/hash.h>
37 #include <libpspp/magic.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/group.h>
44 #include <math/levene.h>
45 #include <output/manager.h>
46 #include <output/table.h>
47 #include "sort-criteria.h"
48
49 #include "gettext.h"
50 #define _(msgid) gettext (msgid)
51
52 /* (headers) */
53
54 /* (specification)
55    "ONEWAY" (oneway_):
56    *^variables=custom;
57    missing=miss:!analysis/listwise,
58            incl:include/!exclude;
59    +contrast= double list;
60    +statistics[st_]=descriptives,homogeneity.
61 */
62 /* (declarations) */
63 /* (functions) */
64
65 static struct cmd_oneway cmd;
66
67 /* The independent variable */
68 static const struct variable *indep_var;
69
70 /* Number of dependent variables */
71 static size_t n_vars;
72
73 /* The dependent variables */
74 static const struct variable **vars;
75
76
77 /* A  hash table containing all the distinct values of the independent
78    variables */
79 static struct hsh_table *global_group_hash ;
80
81 /* The number of distinct values of the independent variable, when all
82    missing values are disregarded */
83 static int ostensible_number_of_groups = -1;
84
85
86 static void run_oneway (struct cmd_oneway *, struct casereader *,
87                         const struct dataset *);
88
89
90 /* Routines to show the output tables */
91 static void show_anova_table(void);
92 static void show_descriptives(void);
93 static void show_homogeneity(void);
94
95 static void show_contrast_coeffs(short *);
96 static void show_contrast_tests(short *);
97
98
99 enum stat_table_t {STAT_DESC = 1, STAT_HOMO = 2};
100
101 static enum stat_table_t stat_tables ;
102
103 void output_oneway(void);
104
105
106 int
107 cmd_oneway (struct lexer *lexer, struct dataset *ds)
108 {
109   struct casegrouper *grouper;
110   struct casereader *group;
111   int i;
112   bool ok;
113
114   if ( !parse_oneway (lexer, ds, &cmd, NULL) )
115     return CMD_FAILURE;
116
117   /* What statistics were requested */
118   if ( cmd.sbc_statistics )
119     {
120
121       for (i = 0 ; i < ONEWAY_ST_count ; ++i )
122         {
123           if  ( ! cmd.a_statistics[i]  ) continue;
124
125           switch (i) {
126           case ONEWAY_ST_DESCRIPTIVES:
127             stat_tables |= STAT_DESC;
128             break;
129           case ONEWAY_ST_HOMOGENEITY:
130             stat_tables |= STAT_HOMO;
131             break;
132           }
133         }
134     }
135
136   /* Data pass.  FIXME: error handling. */
137   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
138   while (casegrouper_get_next_group (grouper, &group))
139     run_oneway (&cmd, group, ds);
140   ok = casegrouper_destroy (grouper);
141   ok = proc_commit (ds) && ok;
142
143   free (vars);
144   free_oneway (&cmd);
145
146   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
147 }
148
149
150 void
151 output_oneway(void)
152 {
153   size_t i;
154   short *bad_contrast ;
155
156   bad_contrast = xnmalloc (cmd.sbc_contrast, sizeof *bad_contrast);
157
158   /* Check the sanity of the given contrast values */
159   for (i = 0 ; i < cmd.sbc_contrast ; ++i )
160     {
161       int j;
162       double sum = 0;
163
164       bad_contrast[i] = 0;
165       if ( subc_list_double_count(&cmd.dl_contrast[i]) !=
166            ostensible_number_of_groups )
167         {
168           msg(SW,
169               _("Number of contrast coefficients must equal the number of groups"));
170           bad_contrast[i] = 1;
171           continue;
172         }
173
174       for (j=0; j < ostensible_number_of_groups ; ++j )
175         sum += subc_list_double_at(&cmd.dl_contrast[i],j);
176
177       if ( sum != 0.0 )
178         msg(SW,_("Coefficients for contrast %d do not total zero"),
179             (int) i + 1);
180     }
181
182   if ( stat_tables & STAT_DESC )
183     show_descriptives();
184
185   if ( stat_tables & STAT_HOMO )
186     show_homogeneity();
187
188   show_anova_table();
189
190   if (cmd.sbc_contrast )
191     {
192       show_contrast_coeffs(bad_contrast);
193       show_contrast_tests(bad_contrast);
194     }
195
196
197   free(bad_contrast);
198
199   /* Clean up */
200   for (i = 0 ; i < n_vars ; ++i )
201     {
202       struct hsh_table *group_hash = group_proc_get (vars[i])->group_hash;
203
204       hsh_destroy(group_hash);
205     }
206
207   hsh_destroy(global_group_hash);
208
209 }
210
211
212
213
214 /* Parser for the variables sub command */
215 static int
216 oneway_custom_variables (struct lexer *lexer,
217                         struct dataset *ds, struct cmd_oneway *cmd UNUSED,
218                         void *aux UNUSED)
219 {
220   struct dictionary *dict = dataset_dict (ds);
221
222   lex_match (lexer, '=');
223
224   if ((lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
225       && lex_token (lexer) != T_ALL)
226     return 2;
227
228   if (!parse_variables_const (lexer, dict, &vars, &n_vars,
229                         PV_DUPLICATE
230                         | PV_NUMERIC | PV_NO_SCRATCH) )
231     {
232       free (vars);
233       return 0;
234     }
235
236   assert(n_vars);
237
238   if ( ! lex_match (lexer, T_BY))
239     return 2;
240
241   indep_var = parse_variable (lexer, dict);
242
243   if ( !indep_var )
244     {
245       msg(SE,_("`%s' is not a variable name"),lex_tokid (lexer));
246       return 0;
247     }
248
249   return 1;
250 }
251
252
253 /* Show the ANOVA table */
254 static void
255 show_anova_table(void)
256 {
257   size_t i;
258   int n_cols =7;
259   size_t n_rows = n_vars * 3 + 1;
260
261   struct tab_table *t;
262
263
264   t = tab_create (n_cols,n_rows,0);
265   tab_headers (t, 2, 0, 1, 0);
266   tab_dim (t, tab_natural_dimensions);
267
268
269   tab_box (t,
270            TAL_2, TAL_2,
271            -1, TAL_1,
272            0, 0,
273            n_cols - 1, n_rows - 1);
274
275   tab_hline (t, TAL_2, 0, n_cols - 1, 1 );
276   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
277   tab_vline (t, TAL_0, 1, 0, 0);
278
279   tab_text (t, 2, 0, TAB_CENTER | TAT_TITLE, _("Sum of Squares"));
280   tab_text (t, 3, 0, TAB_CENTER | TAT_TITLE, _("df"));
281   tab_text (t, 4, 0, TAB_CENTER | TAT_TITLE, _("Mean Square"));
282   tab_text (t, 5, 0, TAB_CENTER | TAT_TITLE, _("F"));
283   tab_text (t, 6, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
284
285
286   for ( i=0 ; i < n_vars ; ++i )
287     {
288       struct group_statistics *totals = &group_proc_get (vars[i])->ugs;
289       struct hsh_table *group_hash = group_proc_get (vars[i])->group_hash;
290       struct hsh_iterator g;
291       struct group_statistics *gs;
292       double ssa=0;
293       const char *s = var_to_string(vars[i]);
294
295       for (gs =  hsh_first (group_hash,&g);
296            gs != 0;
297            gs = hsh_next(group_hash,&g))
298         {
299           ssa += (gs->sum * gs->sum)/gs->n;
300         }
301
302       ssa -= ( totals->sum * totals->sum ) / totals->n ;
303
304       tab_text (t, 0, i * 3 + 1, TAB_LEFT | TAT_TITLE, s);
305       tab_text (t, 1, i * 3 + 1, TAB_LEFT | TAT_TITLE, _("Between Groups"));
306       tab_text (t, 1, i * 3 + 2, TAB_LEFT | TAT_TITLE, _("Within Groups"));
307       tab_text (t, 1, i * 3 + 3, TAB_LEFT | TAT_TITLE, _("Total"));
308
309       if (i > 0)
310         tab_hline(t, TAL_1, 0, n_cols - 1 , i * 3 + 1);
311
312       {
313         struct group_proc *gp = group_proc_get (vars[i]);
314         const double sst = totals->ssq - ( totals->sum * totals->sum) / totals->n ;
315         const double df1 = gp->n_groups - 1;
316         const double df2 = totals->n - gp->n_groups ;
317         const double msa = ssa / df1;
318
319         gp->mse  = (sst - ssa) / df2;
320
321
322         /* Sums of Squares */
323         tab_float (t, 2, i * 3 + 1, 0, ssa, 10, 2);
324         tab_float (t, 2, i * 3 + 3, 0, sst, 10, 2);
325         tab_float (t, 2, i * 3 + 2, 0, sst - ssa, 10, 2);
326
327
328         /* Degrees of freedom */
329         tab_float (t, 3, i * 3 + 1, 0, df1, 4, 0);
330         tab_float (t, 3, i * 3 + 2, 0, df2, 4, 0);
331         tab_float (t, 3, i * 3 + 3, 0, totals->n - 1, 4, 0);
332
333         /* Mean Squares */
334         tab_float (t, 4, i * 3 + 1, TAB_RIGHT, msa, 8, 3);
335         tab_float (t, 4, i * 3 + 2, TAB_RIGHT, gp->mse, 8, 3);
336
337
338         {
339           const double F = msa/gp->mse ;
340
341           /* The F value */
342           tab_float (t, 5, i * 3 + 1, 0,  F, 8, 3);
343
344           /* The significance */
345           tab_float (t, 6, i * 3 + 1, 0, gsl_cdf_fdist_Q(F,df1,df2), 8, 3);
346         }
347
348       }
349
350     }
351
352
353   tab_title (t, _("ANOVA"));
354   tab_submit (t);
355 }
356
357
358 /* Show the descriptives table */
359 static void
360 show_descriptives(void)
361 {
362   size_t v;
363   int n_cols =10;
364   struct tab_table *t;
365   int row;
366
367   const double confidence=0.95;
368   const double q = (1.0 - confidence) / 2.0;
369
370
371   int n_rows = 2 ;
372
373   for ( v = 0 ; v < n_vars ; ++v )
374     n_rows += group_proc_get (vars[v])->n_groups + 1;
375
376   t = tab_create (n_cols,n_rows,0);
377   tab_headers (t, 2, 0, 2, 0);
378   tab_dim (t, tab_natural_dimensions);
379
380
381   /* Put a frame around the entire box, and vertical lines inside */
382   tab_box (t,
383            TAL_2, TAL_2,
384            -1, TAL_1,
385            0, 0,
386            n_cols - 1, n_rows - 1);
387
388   /* Underline headers */
389   tab_hline (t, TAL_2, 0, n_cols - 1, 2 );
390   tab_vline (t, TAL_2, 2, 0, n_rows - 1);
391
392   tab_text (t, 2, 1, TAB_CENTER | TAT_TITLE, _("N"));
393   tab_text (t, 3, 1, TAB_CENTER | TAT_TITLE, _("Mean"));
394   tab_text (t, 4, 1, TAB_CENTER | TAT_TITLE, _("Std. Deviation"));
395   tab_text (t, 5, 1, TAB_CENTER | TAT_TITLE, _("Std. Error"));
396
397
398   tab_vline(t, TAL_0, 7, 0, 0);
399   tab_hline(t, TAL_1, 6, 7, 1);
400   tab_joint_text (t, 6, 0, 7, 0, TAB_CENTER | TAT_TITLE | TAT_PRINTF, _("%g%% Confidence Interval for Mean"),confidence*100.0);
401
402   tab_text (t, 6, 1, TAB_CENTER | TAT_TITLE, _("Lower Bound"));
403   tab_text (t, 7, 1, TAB_CENTER | TAT_TITLE, _("Upper Bound"));
404
405   tab_text (t, 8, 1, TAB_CENTER | TAT_TITLE, _("Minimum"));
406   tab_text (t, 9, 1, TAB_CENTER | TAT_TITLE, _("Maximum"));
407
408
409   tab_title (t, _("Descriptives"));
410
411
412   row = 2;
413   for ( v=0 ; v < n_vars ; ++v )
414     {
415       double T;
416       double std_error;
417
418       struct group_proc *gp = group_proc_get (vars[v]);
419
420       struct group_statistics *gs;
421       struct group_statistics *totals = &gp->ugs;
422
423       const char *s = var_to_string(vars[v]);
424
425       struct group_statistics *const *gs_array =
426         (struct group_statistics *const *) hsh_sort(gp->group_hash);
427       int count = 0;
428
429       tab_text (t, 0, row, TAB_LEFT | TAT_TITLE, s);
430       if ( v > 0)
431         tab_hline(t, TAL_1, 0, n_cols - 1 , row);
432
433       for (count = 0 ; count < hsh_count(gp->group_hash) ; ++count)
434         {
435           gs = gs_array[count];
436
437           tab_text (t, 1, row + count,
438                     TAB_LEFT | TAT_TITLE, var_get_value_name(indep_var,
439                                                              &gs->id));
440
441           /* Now fill in the numbers ... */
442
443           tab_float (t, 2, row + count, 0, gs->n, 8,0);
444
445           tab_float (t, 3, row + count, 0, gs->mean,8,2);
446
447           tab_float (t, 4, row + count, 0, gs->std_dev,8,2);
448
449           std_error = gs->std_dev/sqrt(gs->n) ;
450           tab_float (t, 5, row + count, 0,
451                      std_error, 8,2);
452
453           /* Now the confidence interval */
454
455           T = gsl_cdf_tdist_Qinv(q,gs->n - 1);
456
457           tab_float(t, 6, row + count, 0,
458                     gs->mean - T * std_error, 8, 2);
459
460           tab_float(t, 7, row + count, 0,
461                     gs->mean + T * std_error, 8, 2);
462
463           /* Min and Max */
464
465           tab_float(t, 8, row + count, 0,  gs->minimum, 8, 2);
466           tab_float(t, 9, row + count, 0,  gs->maximum, 8, 2);
467
468         }
469
470       tab_text (t, 1, row + count,
471                 TAB_LEFT | TAT_TITLE ,_("Total"));
472
473       tab_float (t, 2, row + count, 0, totals->n, 8,0);
474
475       tab_float (t, 3, row + count, 0, totals->mean, 8,2);
476
477       tab_float (t, 4, row + count, 0, totals->std_dev,8,2);
478
479       std_error = totals->std_dev/sqrt(totals->n) ;
480
481       tab_float (t, 5, row + count, 0, std_error, 8,2);
482
483       /* Now the confidence interval */
484
485       T = gsl_cdf_tdist_Qinv(q,totals->n - 1);
486
487       tab_float(t, 6, row + count, 0,
488                 totals->mean - T * std_error, 8, 2);
489
490       tab_float(t, 7, row + count, 0,
491                 totals->mean + T * std_error, 8, 2);
492
493       /* Min and Max */
494
495       tab_float(t, 8, row + count, 0,  totals->minimum, 8, 2);
496       tab_float(t, 9, row + count, 0,  totals->maximum, 8, 2);
497
498       row += gp->n_groups + 1;
499     }
500
501
502   tab_submit (t);
503
504
505 }
506
507 /* Show the homogeneity table */
508 static void
509 show_homogeneity(void)
510 {
511   size_t v;
512   int n_cols = 5;
513   size_t n_rows = n_vars + 1;
514
515   struct tab_table *t;
516
517
518   t = tab_create (n_cols,n_rows,0);
519   tab_headers (t, 1, 0, 1, 0);
520   tab_dim (t, tab_natural_dimensions);
521
522   /* Put a frame around the entire box, and vertical lines inside */
523   tab_box (t,
524            TAL_2, TAL_2,
525            -1, TAL_1,
526            0, 0,
527            n_cols - 1, n_rows - 1);
528
529
530   tab_hline(t, TAL_2, 0, n_cols - 1, 1);
531   tab_vline(t, TAL_2, 1, 0, n_rows - 1);
532
533
534   tab_text (t,  1, 0, TAB_CENTER | TAT_TITLE, _("Levene Statistic"));
535   tab_text (t,  2, 0, TAB_CENTER | TAT_TITLE, _("df1"));
536   tab_text (t,  3, 0, TAB_CENTER | TAT_TITLE, _("df2"));
537   tab_text (t,  4, 0, TAB_CENTER | TAT_TITLE, _("Significance"));
538
539
540   tab_title (t, _("Test of Homogeneity of Variances"));
541
542   for ( v=0 ; v < n_vars ; ++v )
543     {
544       double F;
545       const struct variable *var = vars[v];
546       const struct group_proc *gp = group_proc_get (vars[v]);
547       const char *s = var_to_string(var);
548       const struct group_statistics *totals = &gp->ugs;
549
550       const double df1 = gp->n_groups - 1;
551       const double df2 = totals->n - gp->n_groups ;
552
553       tab_text (t, 0, v + 1, TAB_LEFT | TAT_TITLE, s);
554
555       F = gp->levene;
556       tab_float (t, 1, v + 1, TAB_RIGHT, F, 8,3);
557       tab_float (t, 2, v + 1, TAB_RIGHT, df1 ,8,0);
558       tab_float (t, 3, v + 1, TAB_RIGHT, df2 ,8,0);
559
560       /* Now the significance */
561       tab_float (t, 4, v + 1, TAB_RIGHT,gsl_cdf_fdist_Q(F,df1,df2), 8, 3);
562     }
563
564   tab_submit (t);
565 }
566
567
568 /* Show the contrast coefficients table */
569 static void
570 show_contrast_coeffs (short *bad_contrast)
571 {
572   int n_cols = 2 + ostensible_number_of_groups;
573   int n_rows = 2 + cmd.sbc_contrast;
574   union value *group_value;
575   int count = 0 ;
576   void *const *group_values ;
577
578   struct tab_table *t;
579
580   t = tab_create (n_cols,n_rows,0);
581   tab_headers (t, 2, 0, 2, 0);
582   tab_dim (t, tab_natural_dimensions);
583
584   /* Put a frame around the entire box, and vertical lines inside */
585   tab_box (t,
586            TAL_2, TAL_2,
587            -1, TAL_1,
588            0, 0,
589            n_cols - 1, n_rows - 1);
590
591   tab_box (t,
592            -1,-1,
593            TAL_0, TAL_0,
594            2, 0,
595            n_cols - 1, 0);
596
597   tab_box (t,
598            -1,-1,
599            TAL_0, TAL_0,
600            0,0,
601            1,1);
602
603   tab_hline(t, TAL_1, 2, n_cols - 1, 1);
604   tab_hline(t, TAL_2, 0, n_cols - 1, 2);
605
606   tab_vline(t, TAL_2, 2, 0, n_rows - 1);
607
608   tab_title (t, _("Contrast Coefficients"));
609
610   tab_text (t,  0, 2, TAB_LEFT | TAT_TITLE, _("Contrast"));
611
612
613   tab_joint_text (t, 2, 0, n_cols - 1, 0, TAB_CENTER | TAT_TITLE,
614                   var_to_string(indep_var));
615
616   group_values = hsh_sort(global_group_hash);
617   for (count = 0 ;
618        count < hsh_count(global_group_hash) ;
619        ++count)
620     {
621       int i;
622       group_value = group_values[count];
623
624       tab_text (t, count + 2, 1, TAB_CENTER | TAT_TITLE,
625                 var_get_value_name (indep_var, group_value));
626
627       for (i = 0 ; i < cmd.sbc_contrast ; ++i )
628         {
629           tab_text(t, 1, i + 2, TAB_CENTER | TAT_PRINTF, "%d", i + 1);
630
631           if ( bad_contrast[i] )
632             tab_text(t, count + 2, i + 2, TAB_RIGHT, "?" );
633           else
634             tab_text(t, count + 2, i + 2, TAB_RIGHT | TAT_PRINTF, "%g",
635                      subc_list_double_at(&cmd.dl_contrast[i], count)
636                      );
637         }
638     }
639
640   tab_submit (t);
641 }
642
643
644 /* Show the results of the contrast tests */
645 static void
646 show_contrast_tests(short *bad_contrast)
647 {
648   size_t v;
649   int n_cols = 8;
650   size_t n_rows = 1 + n_vars * 2 * cmd.sbc_contrast;
651
652   struct tab_table *t;
653
654   t = tab_create (n_cols,n_rows,0);
655   tab_headers (t, 3, 0, 1, 0);
656   tab_dim (t, tab_natural_dimensions);
657
658   /* Put a frame around the entire box, and vertical lines inside */
659   tab_box (t,
660            TAL_2, TAL_2,
661            -1, TAL_1,
662            0, 0,
663            n_cols - 1, n_rows - 1);
664
665   tab_box (t,
666            -1,-1,
667            TAL_0, TAL_0,
668            0, 0,
669            2, 0);
670
671   tab_hline(t, TAL_2, 0, n_cols - 1, 1);
672   tab_vline(t, TAL_2, 3, 0, n_rows - 1);
673
674
675   tab_title (t, _("Contrast Tests"));
676
677   tab_text (t,  2, 0, TAB_CENTER | TAT_TITLE, _("Contrast"));
678   tab_text (t,  3, 0, TAB_CENTER | TAT_TITLE, _("Value of Contrast"));
679   tab_text (t,  4, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
680   tab_text (t,  5, 0, TAB_CENTER | TAT_TITLE, _("t"));
681   tab_text (t,  6, 0, TAB_CENTER | TAT_TITLE, _("df"));
682   tab_text (t,  7, 0, TAB_CENTER | TAT_TITLE, _("Sig. (2-tailed)"));
683
684   for ( v = 0 ; v < n_vars ; ++v )
685     {
686       int i;
687       int lines_per_variable = 2 * cmd.sbc_contrast;
688
689
690       tab_text (t,  0, (v * lines_per_variable) + 1, TAB_LEFT | TAT_TITLE,
691                 var_to_string(vars[v]));
692
693       for ( i = 0 ; i < cmd.sbc_contrast ; ++i )
694         {
695           int ci;
696           double contrast_value = 0.0;
697           double coef_msq = 0.0;
698           struct group_proc *grp_data = group_proc_get (vars[v]);
699           struct hsh_table *group_hash = grp_data->group_hash;
700
701           void *const *group_stat_array;
702
703           double T;
704           double std_error_contrast ;
705           double df;
706           double sec_vneq=0.0;
707
708
709           /* Note: The calculation of the degrees of freedom in the
710              "variances not equal" case is painfull!!
711              The following formula may help to understand it:
712              \frac{\left(\sum_{i=1}^k{c_i^2\frac{s_i^2}{n_i}}\right)^2}
713              {
714              \sum_{i=1}^k\left(
715              \frac{\left(c_i^2\frac{s_i^2}{n_i}\right)^2}  {n_i-1}
716              \right)
717              }
718           */
719
720           double df_denominator = 0.0;
721           double df_numerator = 0.0;
722           if ( i == 0 )
723             {
724               tab_text (t,  1, (v * lines_per_variable) + i + 1,
725                         TAB_LEFT | TAT_TITLE,
726                         _("Assume equal variances"));
727
728               tab_text (t,  1, (v * lines_per_variable) + i + 1 + cmd.sbc_contrast,
729                         TAB_LEFT | TAT_TITLE,
730                         _("Does not assume equal"));
731             }
732
733           tab_text (t,  2, (v * lines_per_variable) + i + 1,
734                     TAB_CENTER | TAT_TITLE | TAT_PRINTF, "%d",i+1);
735
736
737           tab_text (t,  2, (v * lines_per_variable) + i + 1 + cmd.sbc_contrast,
738                     TAB_CENTER | TAT_TITLE | TAT_PRINTF, "%d",i+1);
739
740
741           if ( bad_contrast[i])
742             continue;
743
744           group_stat_array = hsh_sort(group_hash);
745
746           for (ci = 0 ; ci < hsh_count(group_hash) ;  ++ci)
747             {
748               const double coef = subc_list_double_at(&cmd.dl_contrast[i], ci);
749               struct group_statistics *gs = group_stat_array[ci];
750
751               const double winv = (gs->std_dev * gs->std_dev) / gs->n;
752
753               contrast_value += coef * gs->mean;
754
755               coef_msq += (coef * coef) / gs->n ;
756
757               sec_vneq += (coef * coef) * (gs->std_dev * gs->std_dev ) /gs->n ;
758
759               df_numerator += (coef * coef) * winv;
760               df_denominator += pow2((coef * coef) * winv) / (gs->n - 1);
761             }
762           sec_vneq = sqrt(sec_vneq);
763
764           df_numerator = pow2(df_numerator);
765
766           tab_float (t,  3, (v * lines_per_variable) + i + 1,
767                      TAB_RIGHT, contrast_value, 8,2);
768
769           tab_float (t,  3, (v * lines_per_variable) + i + 1 +
770                      cmd.sbc_contrast,
771                      TAB_RIGHT, contrast_value, 8,2);
772
773           std_error_contrast = sqrt(grp_data->mse * coef_msq);
774
775           /* Std. Error */
776           tab_float (t,  4, (v * lines_per_variable) + i + 1,
777                      TAB_RIGHT, std_error_contrast,
778                      8,3);
779
780           T = fabs(contrast_value / std_error_contrast) ;
781
782           /* T Statistic */
783
784           tab_float (t,  5, (v * lines_per_variable) + i + 1,
785                      TAB_RIGHT, T,
786                      8,3);
787
788           df = grp_data->ugs.n - grp_data->n_groups;
789
790           /* Degrees of Freedom */
791           tab_float (t,  6, (v * lines_per_variable) + i + 1,
792                      TAB_RIGHT,  df,
793                      8,0);
794
795
796           /* Significance TWO TAILED !!*/
797           tab_float (t,  7, (v * lines_per_variable) + i + 1,
798                      TAB_RIGHT,  2 * gsl_cdf_tdist_Q(T,df),
799                      8,3);
800
801
802           /* Now for the Variances NOT Equal case */
803
804           /* Std. Error */
805           tab_float (t,  4,
806                      (v * lines_per_variable) + i + 1 + cmd.sbc_contrast,
807                      TAB_RIGHT, sec_vneq,
808                      8,3);
809
810
811           T = contrast_value / sec_vneq;
812           tab_float (t,  5,
813                      (v * lines_per_variable) + i + 1 + cmd.sbc_contrast,
814                      TAB_RIGHT, T,
815                      8,3);
816
817
818           df = df_numerator / df_denominator;
819
820           tab_float (t,  6,
821                      (v * lines_per_variable) + i + 1 + cmd.sbc_contrast,
822                      TAB_RIGHT, df,
823                      8,3);
824
825           /* The Significance */
826
827           tab_float (t, 7, (v * lines_per_variable) + i + 1 + cmd.sbc_contrast,
828                      TAB_RIGHT,  2 * gsl_cdf_tdist_Q(T,df),
829                      8,3);
830
831
832         }
833
834       if ( v > 0 )
835         tab_hline(t, TAL_1, 0, n_cols - 1, (v * lines_per_variable) + 1);
836     }
837
838   tab_submit (t);
839
840 }
841
842
843 /* ONEWAY ANOVA Calculations */
844
845 static void  postcalc (  struct cmd_oneway *cmd UNUSED );
846
847 static void  precalc ( struct cmd_oneway *cmd UNUSED );
848
849
850
851 /* Pre calculations */
852 static void
853 precalc ( struct cmd_oneway *cmd UNUSED )
854 {
855   size_t i=0;
856
857   for(i=0; i< n_vars ; ++i)
858     {
859       struct group_proc *gp = group_proc_get (vars[i]);
860       struct group_statistics *totals = &gp->ugs;
861
862       /* Create a hash for each of the dependent variables.
863          The hash contains a group_statistics structure,
864          and is keyed by value of the independent variable */
865
866       gp->group_hash =
867         hsh_create(4,
868                    (hsh_compare_func *) compare_group,
869                    (hsh_hash_func *) hash_group,
870                    (hsh_free_func *) free_group,
871                    (void *) var_get_width (indep_var) );
872
873
874       totals->sum=0;
875       totals->n=0;
876       totals->ssq=0;
877       totals->sum_diff=0;
878       totals->maximum = - DBL_MAX;
879       totals->minimum = DBL_MAX;
880     }
881 }
882
883 static void
884 free_value (void *value_, const void *aux UNUSED)
885 {
886   union value *value = value_;
887   free (value);
888 }
889
890 static void
891 run_oneway (struct cmd_oneway *cmd,
892             struct casereader *input,
893             const struct dataset *ds)
894 {
895   struct taint *taint;
896   struct dictionary *dict = dataset_dict (ds);
897   enum mv_class exclude;
898   struct casereader *reader;
899   struct ccase c;
900
901   if (!casereader_peek (input, 0, &c))
902     {
903       casereader_destroy (input);
904       return;
905     }
906   output_split_file_values (ds, &c);
907   case_destroy (&c);
908
909   taint = taint_clone (casereader_get_taint (input));
910
911   global_group_hash = hsh_create(4,
912                                  (hsh_compare_func *) compare_values,
913                                  (hsh_hash_func *) hash_value,
914                                  free_value,
915                                  (void *) var_get_width (indep_var) );
916
917   precalc(cmd);
918
919   exclude = cmd->incl != ONEWAY_INCLUDE ? MV_ANY : MV_SYSTEM;
920   input = casereader_create_filter_missing (input, &indep_var, 1,
921                                             exclude, NULL);
922   if (cmd->miss == ONEWAY_LISTWISE)
923     input = casereader_create_filter_missing (input, vars, n_vars,
924                                               exclude, NULL);
925   input = casereader_create_filter_weight (input, dict, NULL, NULL);
926
927   reader = casereader_clone (input);
928   for (; casereader_read (reader, &c); case_destroy (&c))
929     {
930       size_t i;
931
932       const double weight = dict_get_case_weight (dict, &c, NULL);
933
934       const union value *indep_val = case_data (&c, indep_var);
935       void **p = hsh_probe (global_group_hash, indep_val);
936       if (*p == NULL)
937         *p = value_dup (indep_val, var_get_width (indep_var));
938
939       for ( i = 0 ; i < n_vars ; ++i )
940         {
941           const struct variable *v = vars[i];
942
943           const union value *val = case_data (&c, v);
944
945           struct group_proc *gp = group_proc_get (vars[i]);
946           struct hsh_table *group_hash = gp->group_hash;
947
948           struct group_statistics *gs;
949
950           gs = hsh_find(group_hash, (void *) indep_val );
951
952           if ( ! gs )
953             {
954               gs = xmalloc (sizeof *gs);
955               gs->id = *indep_val;
956               gs->sum=0;
957               gs->n=0;
958               gs->ssq=0;
959               gs->sum_diff=0;
960               gs->minimum = DBL_MAX;
961               gs->maximum = -DBL_MAX;
962
963               hsh_insert ( group_hash, (void *) gs );
964             }
965
966           if (!var_is_value_missing (v, val, exclude))
967             {
968               struct group_statistics *totals = &gp->ugs;
969
970               totals->n+=weight;
971               totals->sum+=weight * val->f;
972               totals->ssq+=weight * val->f * val->f;
973
974               if ( val->f * weight  < totals->minimum )
975                 totals->minimum = val->f * weight;
976
977               if ( val->f * weight  > totals->maximum )
978                 totals->maximum = val->f * weight;
979
980               gs->n+=weight;
981               gs->sum+=weight * val->f;
982               gs->ssq+=weight * val->f * val->f;
983
984               if ( val->f * weight  < gs->minimum )
985                 gs->minimum = val->f * weight;
986
987               if ( val->f * weight  > gs->maximum )
988                 gs->maximum = val->f * weight;
989             }
990
991           gp->n_groups = hsh_count ( group_hash );
992         }
993
994     }
995   casereader_destroy (reader);
996
997   postcalc(cmd);
998
999
1000   if ( stat_tables & STAT_HOMO )
1001     levene (dict, casereader_clone (input), indep_var, n_vars, vars, exclude);
1002
1003   casereader_destroy (input);
1004
1005   ostensible_number_of_groups = hsh_count (global_group_hash);
1006
1007   if (!taint_has_tainted_successor (taint))
1008     output_oneway();
1009   taint_destroy (taint);
1010 }
1011
1012
1013 /* Post calculations for the ONEWAY command */
1014 void
1015 postcalc (  struct cmd_oneway *cmd UNUSED )
1016 {
1017   size_t i=0;
1018
1019
1020   for(i = 0; i < n_vars ; ++i)
1021     {
1022       struct group_proc *gp = group_proc_get (vars[i]);
1023       struct hsh_table *group_hash = gp->group_hash;
1024       struct group_statistics *totals = &gp->ugs;
1025
1026       struct hsh_iterator g;
1027       struct group_statistics *gs;
1028
1029       for (gs =  hsh_first (group_hash,&g);
1030            gs != 0;
1031            gs = hsh_next(group_hash,&g))
1032         {
1033           gs->mean=gs->sum / gs->n;
1034           gs->s_std_dev= sqrt(
1035                               ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1036                               ) ;
1037
1038           gs->std_dev= sqrt(
1039                             gs->n/(gs->n-1) *
1040                             ( (gs->ssq / gs->n ) - gs->mean * gs->mean )
1041                             ) ;
1042
1043           gs->se_mean = gs->std_dev / sqrt(gs->n);
1044           gs->mean_diff= gs->sum_diff / gs->n;
1045
1046         }
1047
1048
1049
1050       totals->mean = totals->sum / totals->n;
1051       totals->std_dev= sqrt(
1052                             totals->n/(totals->n-1) *
1053                             ( (totals->ssq / totals->n ) - totals->mean * totals->mean )
1054                             ) ;
1055
1056       totals->se_mean = totals->std_dev / sqrt(totals->n);
1057
1058     }
1059 }
1060
1061 /*
1062   Local Variables:
1063   mode: c
1064   End:
1065 */