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