f8d8ad36ea8821c281e8cf9d6e38a1117309d3d6
[pspp] / src / examine.q
1 /* PSPP - EXAMINE data for normality . -*-c-*-
2
3 Copyright (C) 2004 Free Software Foundation, Inc.
4 Author: John Darrington 2004
5
6 This program is free software; you can redistribute it and/or
7 modify it under the terms of the GNU General Public License as
8 published by the Free Software Foundation; either version 2 of the
9 License, or (at your option) any later version.
10
11 This program is distributed in the hope that it will be useful, but
12 WITHOUT ANY WARRANTY; without even the implied warranty of
13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
14 General Public License for more details.
15
16 You should have received a copy of the GNU General Public License
17 along with this program; if not, write to the Free Software
18 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
19 02111-1307, USA. */
20
21 #include <config.h>
22 #include <gsl/gsl_cdf.h>
23 #include "error.h"
24 #include <stdio.h>
25 #include <stdlib.h>
26 #include <math.h>
27 #include "alloc.h"
28 #include "str.h"
29 #include "case.h"
30 #include "dictionary.h"
31 #include "command.h"
32 #include "lexer.h"
33 #include "error.h"
34 #include "magic.h"
35 #include "misc.h"
36 #include "tab.h"
37 #include "som.h"
38 #include "value-labels.h"
39 #include "var.h"
40 #include "vfm.h"
41 #include "hash.h"
42 #include "casefile.h"
43 #include "factor_stats.h"
44 /* (headers) */
45 #include "chart.h"
46
47 /* (specification)
48    "EXAMINE" (xmn_):
49    *variables=custom;
50    +total=custom;
51    +nototal=custom;
52    +missing=miss:pairwise/!listwise,
53    rep:report/!noreport,
54    incl:include/!exclude;
55    +compare=cmp:variables/!groups;
56    +plot[plt_]=stemleaf,boxplot,npplot,:spreadlevel(*d:n),histogram,all,none;
57    +cinterval=double;
58    +statistics[st_]=descriptives,:extreme(*d:n),all,none.
59 */
60
61 /* (declarations) */
62
63 /* (functions) */
64
65
66
67 static struct cmd_examine cmd;
68
69 static struct variable **dependent_vars;
70
71 static int n_dependent_vars;
72
73
74 struct factor 
75 {
76   /* The independent variable */
77   struct variable *indep_var[2];
78
79
80   /* Hash table of factor stats indexed by 2 values */
81   struct hsh_table *fstats;
82
83   /* The hash table after it has been crunched */
84   struct factor_statistics **fs;
85
86   struct factor *next;
87
88 };
89
90 /* Linked list of factors */
91 static struct factor *factors=0;
92
93 static struct metrics *totals=0;
94
95 void
96 print_factors(void)
97 {
98   struct factor *f = factors;
99
100   while (f) 
101     {
102       struct  factor_statistics **fs = f->fs;
103
104       printf("Factor: %s BY %s\n", 
105              var_to_string(f->indep_var[0]),
106              var_to_string(f->indep_var[1]) );
107
108
109       printf("Contains %d entries\n", hsh_count(f->fstats));
110
111       
112       while (*fs) 
113         {
114           printf("Factor %g; %g\n", (*fs)->id[0].f, (*fs)->id[1].f);
115           
116           /* 
117              printf("Factor %s; %s\n",
118              value_to_string(&(*fs)->id[0], f->indep_var[0]),
119              value_to_string(&(*fs)->id[1], f->indep_var[1]));
120           */
121
122                  
123           printf("Sum is %g; ",(*fs)->m[0].sum);
124           printf("N is %g; ",(*fs)->m[0].n);
125           printf("Mean is %g\n",(*fs)->m[0].mean);
126
127           fs++ ;
128         }
129
130       f = f->next;
131     }
132
133   
134 }
135
136
137 /* Parse the clause specifying the factors */
138 static int examine_parse_independent_vars(struct cmd_examine *cmd);
139
140
141
142 /* Output functions */
143 static void show_summary(struct variable **dependent_var, int n_dep_var, 
144                          const struct factor *f);
145
146 static void show_extremes(struct variable **dependent_var, 
147                           int n_dep_var, 
148                           const struct factor *factor,
149                           int n_extremities);
150
151 static void show_descriptives(struct variable **dependent_var, 
152                               int n_dep_var, 
153                               struct factor *factor);
154
155
156 void np_plot(const struct metrics *m, const char *factorname);
157
158
159
160
161 /* Per Split function */
162 static void run_examine(const struct casefile *cf, void *cmd_);
163
164 static void output_examine(void);
165
166
167 void factor_calc(struct ccase *c, int case_no, 
168                  double weight, int case_missing);
169
170
171 /* Function to use for testing for missing values */
172 static is_missing_func value_is_missing;
173
174
175 int
176 cmd_examine(void)
177 {
178
179   if ( !parse_examine(&cmd) )
180     return CMD_FAILURE;
181
182   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
183   if (cmd.incl == XMN_INCLUDE ) 
184     value_is_missing = is_system_missing;
185   else
186     value_is_missing = is_missing;
187
188   if ( cmd.st_n == SYSMIS ) 
189     cmd.st_n = 5;
190
191   if ( ! cmd.sbc_cinterval) 
192     cmd.n_cinterval[0] = 95.0;
193
194   multipass_procedure_with_splits (run_examine, &cmd);
195
196   if ( totals ) 
197     free(totals);
198
199   return CMD_SUCCESS;
200 };
201
202
203
204 /* Show all the appropriate tables */
205 static void
206 output_examine(void)
207 {
208   struct factor *fctr;
209
210   /* Show totals if appropriate */
211   if ( ! cmd.sbc_nototal || factors == 0 )
212     {
213       show_summary(dependent_vars, n_dependent_vars, 0);
214
215       if ( cmd.sbc_statistics ) 
216         {
217           if ( cmd.a_statistics[XMN_ST_EXTREME]) 
218             show_extremes(dependent_vars, n_dependent_vars, 0, cmd.st_n);
219
220           if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES]) 
221             show_descriptives(dependent_vars, n_dependent_vars, 0);
222
223         }
224
225       if ( cmd.sbc_plot) 
226         {
227           if ( cmd.a_plot[XMN_PLT_NPPLOT] ) 
228             {
229               int v;
230
231               for ( v = 0 ; v < n_dependent_vars; ++v ) 
232                   np_plot(&totals[v], var_to_string(dependent_vars[v]));
233             }
234         }
235
236
237     }
238
239
240   /* Show grouped statistics  as appropriate */
241   fctr = factors;
242   while ( fctr ) 
243     {
244       show_summary(dependent_vars, n_dependent_vars, fctr);
245
246       if ( cmd.sbc_statistics ) 
247         {
248           if ( cmd.a_statistics[XMN_ST_EXTREME]) 
249             show_extremes(dependent_vars, n_dependent_vars, fctr, cmd.st_n);
250
251           if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES]) 
252             show_descriptives(dependent_vars, n_dependent_vars, fctr);
253         }
254
255       if ( cmd.sbc_plot) 
256         {
257           if ( cmd.a_plot[XMN_PLT_NPPLOT] ) 
258             {
259               int v;
260               for ( v = 0 ; v < n_dependent_vars; ++ v)
261                 {
262                   
263                   struct factor_statistics **fs = fctr->fs ;
264                   for ( fs = fctr->fs ; *fs ; ++fs ) 
265                     {
266                       char buf1[100];
267                       char buf2[100];
268                       sprintf(buf1, "%s (",
269                               var_to_string(dependent_vars[v]));
270                       
271                       sprintf(buf2, "%s = %s",
272                              var_to_string(fctr->indep_var[0]),
273                              value_to_string(&(*fs)->id[0],fctr->indep_var[0]));
274                       
275                       strcat(buf1, buf2);
276
277                       
278                       if ( fctr->indep_var[1] ) 
279                         {
280                           sprintf(buf2, "; %s = %s)",
281                                   var_to_string(fctr->indep_var[1]),
282                                   value_to_string(&(*fs)->id[1],
283                                                   fctr->indep_var[1]));
284                           strcat(buf1, buf2);
285                         }
286                       else
287                         {
288                           strcat(buf1, ")");
289                         }
290
291                       np_plot(&(*fs)->m[v],buf1);
292
293                     }
294                   
295                 }
296
297             }
298         }
299
300       fctr = fctr->next;
301     }
302
303 }
304
305
306
307 /* TOTAL and NOTOTAL are simple, mutually exclusive flags */
308 static int
309 xmn_custom_total(struct cmd_examine *p)
310 {
311   if ( p->sbc_nototal ) 
312     {
313       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
314       return 0;
315     }
316
317   return 1;
318 }
319
320 static int
321 xmn_custom_nototal(struct cmd_examine *p)
322 {
323   if ( p->sbc_total ) 
324     {
325       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
326       return 0;
327     }
328
329   return 1;
330 }
331
332
333
334 /* Parser for the variables sub command */
335 static int
336 xmn_custom_variables(struct cmd_examine *cmd )
337 {
338
339   lex_match('=');
340
341   if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
342       && token != T_ALL)
343     return 2;
344   
345   if (!parse_variables (default_dict, &dependent_vars, &n_dependent_vars,
346                         PV_NO_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH) )
347     {
348       free (dependent_vars);
349       return 0;
350     }
351
352   assert(n_dependent_vars);
353
354   totals = xmalloc( sizeof(struct metrics) * n_dependent_vars);
355
356   if ( lex_match(T_BY))
357     {
358       return examine_parse_independent_vars(cmd);
359     }
360
361   return 1;
362 }
363
364
365
366 /* Parse the clause specifying the factors */
367 static int
368 examine_parse_independent_vars(struct cmd_examine *cmd)
369 {
370
371   struct factor *sf = xmalloc(sizeof(struct factor));
372
373   if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
374       && token != T_ALL)
375     return 2;
376
377
378   sf->indep_var[0] = parse_variable();
379   sf->indep_var[1] = 0;
380
381   if ( token == T_BY ) 
382     {
383
384       lex_match(T_BY);
385
386       if ((token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
387           && token != T_ALL)
388         return 2;
389
390       sf->indep_var[1] = parse_variable();
391
392     }
393
394
395   sf->fstats = hsh_create(4,
396                           (hsh_compare_func *) factor_statistics_compare,
397                           (hsh_hash_func *) factor_statistics_hash,
398                           (hsh_free_func *) factor_statistics_free,
399                           0);
400
401   sf->next = factors;
402   factors = sf;
403   
404   lex_match(',');
405
406   if ( token == '.' || token == '/' ) 
407     return 1;
408
409   return examine_parse_independent_vars(cmd);
410 }
411
412
413
414
415 void populate_descriptives(struct tab_table *t, int col, int row, 
416                            const struct metrics *fs);
417
418 void populate_extremes(struct tab_table *t, int col, int row, int n, 
419                        const struct metrics *m);
420
421 void populate_summary(struct tab_table *t, int col, int row,
422                       const struct metrics *m);
423
424
425
426
427 static int bad_weight_warn = 1;
428
429
430 /* Perform calculations for the sub factors */
431 void
432 factor_calc(struct ccase *c, int case_no, double weight, int case_missing)
433 {
434   int v;
435   struct factor *fctr = factors;
436
437   while ( fctr) 
438     {
439       union value indep_vals[2] ;
440
441       indep_vals[0] = * case_data(c, fctr->indep_var[0]->fv);
442
443       if ( fctr->indep_var[1] ) 
444         indep_vals[1] = * case_data(c, fctr->indep_var[1]->fv);
445       else
446         indep_vals[1].f = SYSMIS;
447
448       assert(fctr->fstats);
449
450       struct factor_statistics **foo = ( struct factor_statistics ** ) 
451         hsh_probe(fctr->fstats, (void *) &indep_vals);
452
453       if ( !*foo ) 
454         {
455
456           *foo = create_factor_statistics(n_dependent_vars, 
457                                           &indep_vals[0],
458                                           &indep_vals[1]);
459
460           for ( v =  0 ; v  < n_dependent_vars ; ++v ) 
461             {
462               metrics_precalc( &(*foo)->m[v] );
463             }
464
465         }
466
467       for ( v =  0 ; v  < n_dependent_vars ; ++v ) 
468         {
469           const struct variable *var = dependent_vars[v];
470           const union value *val = case_data (c, var->fv);
471
472           if ( value_is_missing(val,var) || case_missing ) 
473             val = 0;
474
475           metrics_calc( &(*foo)->m[v], val, weight, case_no );
476         }
477
478       fctr = fctr->next;
479     }
480
481
482 }
483
484
485
486
487 static void 
488 run_examine(const struct casefile *cf, void *cmd_ )
489 {
490   struct casereader *r;
491   struct ccase c;
492   int v;
493
494   const struct cmd_examine *cmd = (struct cmd_examine *) cmd_;
495
496   /* Make sure we haven't got rubbish left over from a 
497      previous split */
498
499   struct factor *fctr = factors;
500   while (fctr) 
501     {
502       struct factor *next = fctr->next;
503
504       hsh_clear(fctr->fstats);
505
506       fctr->fs = 0;
507
508       fctr = next;
509     }
510
511
512
513   for ( v = 0 ; v < n_dependent_vars ; ++v ) 
514     metrics_precalc(&totals[v]);
515
516   for(r = casefile_get_reader (cf);
517       casereader_read (r, &c) ;
518       case_destroy (&c) ) 
519     {
520       int case_missing=0;
521       const int case_no = casereader_cnum(r);
522
523       const double weight = 
524         dict_get_case_weight(default_dict, &c, &bad_weight_warn);
525
526       if ( cmd->miss == XMN_LISTWISE ) 
527         {
528           for ( v = 0 ; v < n_dependent_vars ; ++v ) 
529             {
530               const struct variable *var = dependent_vars[v];
531               const union value *val = case_data (&c, var->fv);
532
533               if ( value_is_missing(val,var))
534                 case_missing = 1;
535                    
536             }
537         }
538
539       for ( v = 0 ; v < n_dependent_vars ; ++v ) 
540         {
541           const struct variable *var = dependent_vars[v];
542           const union value *val = case_data (&c, var->fv);
543
544           if ( value_is_missing(val,var) || case_missing ) 
545             val = 0;
546
547           metrics_calc(&totals[v], val, weight, case_no );
548     
549         }
550
551       factor_calc(&c, case_no, weight, case_missing);
552
553     }
554
555
556   for ( v = 0 ; v < n_dependent_vars ; ++v)
557     {
558       fctr = factors;
559       while ( fctr ) 
560         {
561           struct hsh_iterator hi;
562           struct factor_statistics *fs;
563
564           for ( fs = hsh_first(fctr->fstats, &hi);
565                 fs != 0 ;
566                 fs = hsh_next(fctr->fstats, &hi))
567             {
568               metrics_postcalc(&fs->m[v]);
569             }
570
571           fctr = fctr->next;
572         }
573       metrics_postcalc(&totals[v]);
574     }
575
576
577   /* Make sure that the combination of factors are complete */
578
579   fctr = factors;
580   while ( fctr ) 
581     {
582       struct hsh_iterator hi;
583       struct hsh_iterator hi0;
584       struct hsh_iterator hi1;
585       struct factor_statistics *fs;
586
587       struct hsh_table *idh0=0;
588       struct hsh_table *idh1=0;
589       union value *val0;
590       union value *val1;
591           
592       idh0 = hsh_create(4, (hsh_compare_func *) compare_values,
593                         (hsh_hash_func *) hash_value,
594                         0,0);
595
596       idh1 = hsh_create(4, (hsh_compare_func *) compare_values,
597                         (hsh_hash_func *) hash_value,
598                         0,0);
599
600
601       for ( fs = hsh_first(fctr->fstats, &hi);
602             fs != 0 ;
603             fs = hsh_next(fctr->fstats, &hi))
604         {
605           hsh_insert(idh0,(void *) &fs->id[0]);
606           hsh_insert(idh1,(void *) &fs->id[1]);
607         }
608
609       /* Ensure that the factors combination is complete */
610       for ( val0 = hsh_first(idh0, &hi0);
611             val0 != 0 ;
612             val0 = hsh_next(idh0, &hi0))
613         {
614           for ( val1 = hsh_first(idh1, &hi1);
615                 val1 != 0 ;
616                 val1 = hsh_next(idh1, &hi1))
617             {
618               struct factor_statistics **ffs;
619               union value key[2];
620               key[0] = *val0;
621               key[1] = *val1;
622                   
623               ffs = (struct factor_statistics **) 
624                 hsh_probe(fctr->fstats, (void *) &key );
625
626               if ( !*ffs ) {
627                 int i;
628                 (*ffs) = create_factor_statistics (n_dependent_vars,
629                                                    &key[0], &key[1]);
630                 for ( i = 0 ; i < n_dependent_vars ; ++i ) 
631                   metrics_precalc( &(*ffs)->m[i]);
632               }
633             }
634         }
635
636       hsh_destroy(idh0);
637       hsh_destroy(idh1);
638
639       fctr->fs = (struct factor_statistics **) hsh_sort_copy(fctr->fstats);
640
641       fctr = fctr->next;
642     }
643
644   output_examine();
645
646   for ( v = 0 ; v < n_dependent_vars ; ++v ) 
647     hsh_destroy(totals[v].ordered_data);
648
649 }
650
651
652 static void
653 show_summary(struct variable **dependent_var, int n_dep_var, 
654              const struct factor *fctr)
655 {
656   static const char *subtitle[]=
657     {
658       N_("Valid"),
659       N_("Missing"),
660       N_("Total")
661     };
662
663   int i;
664   int heading_columns ;
665   int n_cols;
666   const int heading_rows = 3;
667   struct tab_table *tbl;
668
669   int n_rows ;
670   int n_factors = 1;
671
672   if ( fctr )
673     {
674       heading_columns = 2;
675       n_factors = hsh_count(fctr->fstats);
676       n_rows = n_dep_var * n_factors ;
677
678       if ( fctr->indep_var[1] )
679           heading_columns = 3;
680     }
681   else
682     {
683       heading_columns = 1;
684       n_rows = n_dep_var;
685     }
686
687   n_rows += heading_rows;
688
689   n_cols = heading_columns + 6;
690
691   tbl = tab_create (n_cols,n_rows,0);
692   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
693
694   tab_dim (tbl, tab_natural_dimensions);
695   
696   /* Outline the box */
697   tab_box (tbl, 
698            TAL_2, TAL_2,
699            -1, -1,
700            0, 0,
701            n_cols - 1, n_rows - 1);
702
703   /* Vertical lines for the data only */
704   tab_box (tbl, 
705            -1, -1,
706            -1, TAL_1,
707            heading_columns, 0,
708            n_cols - 1, n_rows - 1);
709
710
711   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
712   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, 1 );
713   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, heading_rows -1 );
714
715   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
716
717
718   tab_title (tbl, 0, _("Case Processing Summary"));
719   
720
721   tab_joint_text(tbl, heading_columns, 0, 
722                  n_cols -1, 0,
723                  TAB_CENTER | TAT_TITLE,
724                  _("Cases"));
725
726   /* Remove lines ... */
727   tab_box (tbl, 
728            -1, -1,
729            TAL_0, TAL_0,
730            heading_columns, 0,
731            n_cols - 1, 0);
732
733   for ( i = 0 ; i < 3 ; ++i ) 
734     {
735       tab_text (tbl, heading_columns + i*2 , 2, TAB_CENTER | TAT_TITLE, 
736                 _("N"));
737
738       tab_text (tbl, heading_columns + i*2 + 1, 2, TAB_CENTER | TAT_TITLE, 
739                 _("Percent"));
740
741       tab_joint_text(tbl, heading_columns + i*2 , 1,
742                      heading_columns + i*2 + 1, 1,
743                      TAB_CENTER | TAT_TITLE,
744                      subtitle[i]);
745
746       tab_box (tbl, -1, -1,
747                TAL_0, TAL_0,
748                heading_columns + i*2, 1,
749                heading_columns + i*2 + 1, 1);
750
751     }
752
753
754   /* Titles for the independent variables */
755   if ( fctr ) 
756     {
757       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE, 
758                 var_to_string(fctr->indep_var[0]));
759
760       if ( fctr->indep_var[1] ) 
761         {
762           tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE, 
763                     var_to_string(fctr->indep_var[1]));
764         }
765                 
766     }
767
768
769   for ( i = 0 ; i < n_dep_var ; ++i ) 
770     {
771       int n_factors = 1;
772       if ( fctr ) 
773         n_factors = hsh_count(fctr->fstats);
774       
775
776       if ( i > 0 ) 
777         tab_hline(tbl, TAL_1, 0, n_cols -1 , i * n_factors + heading_rows);
778       
779       tab_text (tbl, 
780                 0, i * n_factors + heading_rows,
781                 TAB_LEFT | TAT_TITLE, 
782                 var_to_string(dependent_var[i])
783                 );
784
785
786       if ( !fctr ) 
787         populate_summary(tbl, heading_columns, 
788                          (i * n_factors) + heading_rows,
789                          &totals[i]);
790
791
792       else
793         {
794           struct factor_statistics **fs = fctr->fs;
795           int count = 0 ;
796
797           while (*fs) 
798             {
799               static union value prev;
800               
801               if ( 0 != compare_values(&prev, &(*fs)->id[0], 
802                                        fctr->indep_var[0]->width))
803                 {
804                    tab_text (tbl, 
805                              1,
806                              (i * n_factors ) + count + 
807                              heading_rows,
808                              TAB_LEFT | TAT_TITLE, 
809                              value_to_string(&(*fs)->id[0], fctr->indep_var[0])
810                              );
811
812                    if (fctr->indep_var[1] && count > 0 ) 
813                      tab_hline(tbl, TAL_1, 1, n_cols - 1, 
814                                (i * n_factors ) + count + heading_rows);
815
816                 }
817               
818               prev = (*fs)->id[0];
819
820
821               if ( fctr->indep_var[1]) 
822                 tab_text (tbl, 
823                           2,
824                           (i * n_factors ) + count + 
825                           heading_rows,
826                           TAB_LEFT | TAT_TITLE, 
827                           value_to_string(&(*fs)->id[1], fctr->indep_var[1])
828                           );
829
830               populate_summary(tbl, heading_columns, 
831                                (i * n_factors) + count 
832                                + heading_rows,
833                                &(*fs)->m[i]);
834
835               count++ ; 
836               fs++;
837             }
838         }
839     }
840
841   tab_submit (tbl);
842 }
843
844
845 void 
846 populate_summary(struct tab_table *t, int col, int row,
847                  const struct metrics *m)
848
849 {
850   const double total = m->n + m->n_missing ; 
851
852   tab_float(t, col + 0, row + 0, TAB_RIGHT, m->n, 8, 0);
853   tab_float(t, col + 2, row + 0, TAB_RIGHT, m->n_missing, 8, 0);
854   tab_float(t, col + 4, row + 0, TAB_RIGHT, total, 8, 0);
855
856
857   if ( total > 0 ) {
858     tab_text (t, col + 1, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%", 
859               100.0 * m->n / total );
860
861     tab_text (t, col + 3, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%", 
862               100.0 * m->n_missing / total );
863
864     /* This seems a bit pointless !!! */
865     tab_text (t, col + 5, row + 0, TAB_RIGHT | TAT_PRINTF, "%2.0f%%", 
866               100.0 * total / total );
867
868
869   }
870
871
872 }  
873
874
875
876 static void 
877 show_extremes(struct variable **dependent_var, int n_dep_var, 
878               const struct factor *fctr, int n_extremities)
879 {
880   int i;
881   int heading_columns ;
882   int n_cols;
883   const int heading_rows = 1;
884   struct tab_table *tbl;
885
886   int n_factors = 1;
887   int n_rows ;
888
889   if ( fctr )
890     {
891       heading_columns = 2;
892       n_factors = hsh_count(fctr->fstats);
893
894       n_rows = n_dep_var * 2 * n_extremities * n_factors;
895
896       if ( fctr->indep_var[1] )
897           heading_columns = 3;
898     }
899   else
900     {
901       heading_columns = 1;
902       n_rows = n_dep_var * 2 * n_extremities;
903     }
904
905   n_rows += heading_rows;
906
907   heading_columns += 2;
908   n_cols = heading_columns + 2;
909
910   tbl = tab_create (n_cols,n_rows,0);
911   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
912
913   tab_dim (tbl, tab_natural_dimensions);
914   
915   /* Outline the box, No internal lines*/
916   tab_box (tbl, 
917            TAL_2, TAL_2,
918            -1, -1,
919            0, 0,
920            n_cols - 1, n_rows - 1);
921
922   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
923
924   tab_title (tbl, 0, _("Extreme Values"));
925
926
927   tab_vline (tbl, TAL_2, n_cols - 2, 0, n_rows -1);
928   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows -1);
929
930   if ( fctr ) 
931     {
932       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE, 
933                 var_to_string(fctr->indep_var[0]));
934
935       if ( fctr->indep_var[1] ) 
936         tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE, 
937                   var_to_string(fctr->indep_var[1]));
938     }
939
940   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _("Value"));
941   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _("Case Number"));
942
943
944
945
946   for ( i = 0 ; i < n_dep_var ; ++i ) 
947     {
948
949       if ( i > 0 ) 
950         tab_hline(tbl, TAL_1, 0, n_cols -1 , 
951                   i * 2 * n_extremities * n_factors + heading_rows);
952       
953       tab_text (tbl, 0,
954                 i * 2 * n_extremities * n_factors  + heading_rows,
955                 TAB_LEFT | TAT_TITLE, 
956                 var_to_string(dependent_var[i])
957                 );
958
959
960       if ( !fctr ) 
961         populate_extremes(tbl, heading_columns - 2, 
962                           i * 2 * n_extremities * n_factors  + heading_rows,
963                           n_extremities, &totals[i]);
964
965       else
966         {
967           struct factor_statistics **fs = fctr->fs;
968           int count = 0 ;
969
970           while (*fs) 
971             {
972               static union value prev ;
973
974               const int row = heading_rows + ( 2 * n_extremities )  * 
975                 ( ( i  * n_factors  ) +  count );
976
977
978               if ( 0 != compare_values(&prev, &(*fs)->id[0], 
979                                        fctr->indep_var[0]->width))
980                 {
981                   
982                   if ( count > 0 ) 
983                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
984
985                   tab_text (tbl, 
986                             1, row,
987                             TAB_LEFT | TAT_TITLE, 
988                             value_to_string(&(*fs)->id[0], fctr->indep_var[0])
989                             );
990                 }
991
992               prev = (*fs)->id[0];
993
994               if (fctr->indep_var[1] && count > 0 ) 
995                 tab_hline(tbl, TAL_1, 2, n_cols - 1, row);
996
997               if ( fctr->indep_var[1]) 
998                 tab_text (tbl, 2, row,
999                           TAB_LEFT | TAT_TITLE, 
1000                           value_to_string(&(*fs)->id[1], fctr->indep_var[1])
1001                           );
1002
1003               populate_extremes(tbl, heading_columns - 2, 
1004                                 row, n_extremities,
1005                                 &(*fs)->m[i]);
1006
1007               count++ ; 
1008               fs++;
1009             }
1010         }
1011     }
1012
1013   tab_submit(tbl);
1014 }
1015
1016
1017
1018 /* Fill in the extremities table */
1019 void 
1020 populate_extremes(struct tab_table *t, 
1021                   int col, int row, int n, const struct metrics *m)
1022 {
1023   int extremity;
1024   int idx=0;
1025
1026
1027   tab_text(t, col, row,
1028            TAB_RIGHT | TAT_TITLE ,
1029            _("Highest")
1030            );
1031
1032   tab_text(t, col, row + n ,
1033            TAB_RIGHT | TAT_TITLE ,
1034            _("Lowest")
1035            );
1036
1037
1038   tab_hline(t, TAL_1, col, col + 3, row + n );
1039             
1040   for (extremity = 0; extremity < n ; ++extremity ) 
1041     {
1042       /* Highest */
1043       tab_float(t, col + 1, row + extremity,
1044                 TAB_RIGHT,
1045                 extremity + 1, 8, 0);
1046
1047
1048       /* Lowest */
1049       tab_float(t, col + 1, row + extremity + n,
1050                 TAB_RIGHT,
1051                 extremity + 1, 8, 0);
1052
1053     }
1054
1055
1056   /* Lowest */
1057   for (idx = 0, extremity = 0; extremity < n && idx < m->n_data ; ++idx ) 
1058     {
1059       int j;
1060       const struct weighted_value *wv = m->wvp[idx];
1061       struct case_node *cn = wv->case_nos;
1062
1063       
1064       for (j = 0 ; j < wv->w ; ++j  )
1065         {
1066           if ( extremity + j >= n ) 
1067             break ;
1068
1069           tab_float(t, col + 3, row + extremity + j  + n,
1070                     TAB_RIGHT,
1071                     wv->v.f, 8, 2);
1072
1073           tab_float(t, col + 2, row + extremity + j  + n,
1074                     TAB_RIGHT,
1075                     cn->num, 8, 0);
1076
1077           if ( cn->next ) 
1078               cn = cn->next;
1079
1080         }
1081
1082       extremity +=  wv->w ;
1083     }
1084
1085
1086   /* Highest */
1087   for (idx = m->n_data - 1, extremity = 0; extremity < n && idx >= 0; --idx ) 
1088     {
1089       int j;
1090       const struct weighted_value *wv = m->wvp[idx];
1091       struct case_node *cn = wv->case_nos;
1092
1093       for (j = 0 ; j < wv->w ; ++j  )
1094         {
1095           if ( extremity + j >= n ) 
1096             break ;
1097
1098           tab_float(t, col + 3, row + extremity + j,
1099                     TAB_RIGHT,
1100                     wv->v.f, 8, 2);
1101
1102           tab_float(t, col + 2, row + extremity + j,
1103                     TAB_RIGHT,
1104                     cn->num, 8, 0);
1105
1106           if ( cn->next ) 
1107               cn = cn->next;
1108
1109         }
1110
1111       extremity +=  wv->w ;
1112     }
1113 }
1114
1115
1116 /* Show the descriptives table */
1117 void
1118 show_descriptives(struct variable **dependent_var, 
1119                   int n_dep_var, 
1120                   struct factor *fctr)
1121 {
1122   int i;
1123   int heading_columns ;
1124   int n_cols;
1125   const int n_stat_rows = 13;
1126
1127   const int heading_rows = 1;
1128
1129   struct tab_table *tbl;
1130
1131   int n_factors = 1;
1132   int n_rows ;
1133
1134   if ( fctr )
1135     {
1136       heading_columns = 4;
1137       n_factors = hsh_count(fctr->fstats);
1138
1139       n_rows = n_dep_var * n_stat_rows * n_factors;
1140
1141       if ( fctr->indep_var[1] )
1142           heading_columns = 5;
1143     }
1144   else
1145     {
1146       heading_columns = 3;
1147       n_rows = n_dep_var * n_stat_rows;
1148     }
1149
1150   n_rows += heading_rows;
1151
1152   n_cols = heading_columns + 2;
1153
1154
1155   tbl = tab_create (n_cols, n_rows, 0);
1156
1157   tab_headers (tbl, heading_columns + 1, 0, heading_rows, 0);
1158
1159   tab_dim (tbl, tab_natural_dimensions);
1160
1161   /* Outline the box and have no internal lines*/
1162   tab_box (tbl, 
1163            TAL_2, TAL_2,
1164            -1, -1,
1165            0, 0,
1166            n_cols - 1, n_rows - 1);
1167
1168   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1169
1170   tab_vline (tbl, TAL_1, heading_columns, 0, n_rows - 1);
1171   tab_vline (tbl, TAL_2, n_cols - 2, 0, n_rows - 1);
1172   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1173
1174   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE, _("Statistic"));
1175   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE, _("Std. Error"));
1176
1177   tab_title (tbl, 0, _("Descriptives"));
1178
1179
1180   for ( i = 0 ; i < n_dep_var ; ++i ) 
1181     {
1182       const int row = heading_rows + i * n_stat_rows * n_factors ;
1183
1184       if ( i > 0 )
1185         tab_hline(tbl, TAL_1, 0, n_cols - 1, row );
1186
1187       tab_text (tbl, 0,
1188                 i * n_stat_rows * n_factors  + heading_rows,
1189                 TAB_LEFT | TAT_TITLE, 
1190                 var_to_string(dependent_var[i])
1191                 );
1192
1193
1194       if ( fctr  )
1195         {
1196           struct factor_statistics **fs = fctr->fs;
1197           int count = 0;
1198
1199           tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE, 
1200                     var_to_string(fctr->indep_var[0]));
1201
1202
1203           if ( fctr->indep_var[1])
1204             tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE, 
1205                       var_to_string(fctr->indep_var[1]));
1206
1207           while( *fs ) 
1208             {
1209
1210               static union value prev ;
1211
1212               const int row = heading_rows + n_stat_rows  * 
1213                 ( ( i  * n_factors  ) +  count );
1214
1215
1216               if ( 0 != compare_values(&prev, &(*fs)->id[0], 
1217                                        fctr->indep_var[0]->width))
1218                 {
1219                   
1220                   if ( count > 0 ) 
1221                     tab_hline (tbl, TAL_1, 1, n_cols - 1, row);
1222
1223                   tab_text (tbl, 
1224                             1, row,
1225                             TAB_LEFT | TAT_TITLE, 
1226                             value_to_string(&(*fs)->id[0], fctr->indep_var[0])
1227                             );
1228                 }
1229
1230               prev = (*fs)->id[0];
1231
1232               if (fctr->indep_var[1] && count > 0 ) 
1233                 tab_hline(tbl, TAL_1, 2, n_cols - 1, row);
1234
1235               if ( fctr->indep_var[1]) 
1236                 tab_text (tbl, 2, row,
1237                           TAB_LEFT | TAT_TITLE, 
1238                           value_to_string(&(*fs)->id[1], fctr->indep_var[1])
1239                           );
1240
1241               populate_descriptives(tbl, heading_columns - 2, 
1242                                 row, &(*fs)->m[i]);
1243
1244               count++ ; 
1245               fs++;
1246             }
1247
1248         }
1249
1250       else 
1251         {
1252           
1253           populate_descriptives(tbl, heading_columns - 2, 
1254                                 i * n_stat_rows * n_factors  + heading_rows,
1255                                 &totals[i]);
1256         }
1257     }
1258
1259   tab_submit(tbl);
1260
1261 }
1262
1263
1264
1265
1266
1267
1268 /* Fill in the descriptives data */
1269 void
1270 populate_descriptives(struct tab_table *tbl, int col, int row, 
1271                       const struct metrics *m)
1272 {
1273
1274   const double t = gsl_cdf_tdist_Qinv(1 - cmd.n_cinterval[0]/100.0/2.0, \
1275                                       m->n -1);
1276
1277
1278   tab_text (tbl, col, 
1279             row,
1280             TAB_LEFT | TAT_TITLE,
1281             _("Mean"));
1282
1283   tab_float (tbl, col + 2,
1284              row,
1285              TAB_CENTER,
1286              m->mean,
1287              8,2);
1288   
1289   tab_float (tbl, col + 3,
1290              row,
1291              TAB_CENTER,
1292              m->stderr,
1293              8,3);
1294   
1295
1296   tab_text (tbl, col, 
1297             row + 1,
1298             TAB_LEFT | TAT_TITLE | TAT_PRINTF,
1299             _("%g%% Confidence Interval for Mean"), cmd.n_cinterval[0]);
1300
1301
1302   tab_text (tbl, col + 1, 
1303             row  + 1,
1304             TAB_LEFT | TAT_TITLE,
1305             _("Lower Bound"));
1306
1307   tab_float (tbl, col + 2,
1308              row + 1,
1309              TAB_CENTER,
1310              m->mean - t * m->stderr, 
1311              8,3);
1312
1313   tab_text (tbl, col + 1,  
1314             row + 2,
1315             TAB_LEFT | TAT_TITLE,
1316             _("Upper Bound"));
1317
1318
1319   tab_float (tbl, col + 2,
1320              row + 2,
1321              TAB_CENTER,
1322              m->mean + t * m->stderr, 
1323              8,3);
1324
1325   tab_text (tbl, col, 
1326             row + 3,
1327             TAB_LEFT | TAT_TITLE,
1328             _("5% Trimmed Mean"));
1329
1330   tab_float (tbl, col + 2, 
1331              row + 3,
1332              TAB_CENTER,
1333              m->trimmed_mean,
1334              8,2);
1335
1336   tab_text (tbl, col, 
1337             row + 4,
1338             TAB_LEFT | TAT_TITLE,
1339             _("Median"));
1340
1341   tab_text (tbl, col, 
1342             row + 5,
1343             TAB_LEFT | TAT_TITLE,
1344             _("Variance"));
1345
1346   tab_float (tbl, col + 2,
1347              row + 5,
1348              TAB_CENTER,
1349              m->var,
1350              8,3);
1351
1352
1353   tab_text (tbl, col, 
1354             row + 6,
1355             TAB_LEFT | TAT_TITLE,
1356             _("Std. Deviation"));
1357
1358
1359   tab_float (tbl, col + 2,
1360              row + 6,
1361              TAB_CENTER,
1362              m->stddev,
1363              8,3);
1364
1365   
1366   tab_text (tbl, col, 
1367             row + 7,
1368             TAB_LEFT | TAT_TITLE,
1369             _("Minimum"));
1370
1371   tab_float (tbl, col + 2,
1372              row + 7,
1373              TAB_CENTER,
1374              m->min,
1375              8,3);
1376
1377   tab_text (tbl, col, 
1378             row + 8,
1379             TAB_LEFT | TAT_TITLE,
1380             _("Maximum"));
1381
1382   tab_float (tbl, col + 2,
1383              row + 8,
1384              TAB_CENTER,
1385              m->max,
1386              8,3);
1387
1388
1389   tab_text (tbl, col, 
1390             row + 9,
1391             TAB_LEFT | TAT_TITLE,
1392             _("Range"));
1393
1394
1395   tab_float (tbl, col + 2,
1396              row + 9,
1397              TAB_CENTER,
1398              m->max - m->min,
1399              8,3);
1400
1401   tab_text (tbl, col, 
1402             row + 10,
1403             TAB_LEFT | TAT_TITLE,
1404             _("Interquartile Range"));
1405
1406   tab_text (tbl, col, 
1407             row + 11,
1408             TAB_LEFT | TAT_TITLE,
1409             _("Skewness"));
1410
1411   tab_text (tbl, col, 
1412             row + 12,
1413             TAB_LEFT | TAT_TITLE,
1414             _("Kurtosis"));
1415 }
1416
1417
1418
1419
1420
1421 /* Plot the normal and detrended normal plots for m
1422    Label the plots with factorname */
1423 void
1424 np_plot(const struct metrics *m, const char *factorname)
1425 {
1426   int i;
1427   double yfirst=0, ylast=0;
1428
1429   /* Normal Plot */
1430   struct chart np_chart;
1431
1432   /* Detrended Normal Plot */
1433   struct chart dnp_chart;
1434
1435   const struct weighted_value *wv = *(m->wvp);
1436
1437
1438   /* The slope and intercept of the ideal normal probability line */
1439   const double slope = 1.0 / m->stddev;
1440   const double intercept = - m->mean / m->stddev;
1441
1442   /* Cowardly refuse to plot an empty data set */
1443   if ( m->n_data == 0 ) 
1444     return ; 
1445
1446   chart_initialise(&np_chart);
1447   chart_write_title(&np_chart, _("Normal Q-Q Plot of %s"), factorname);
1448   chart_write_xlabel(&np_chart, _("Observed Value"));
1449   chart_write_ylabel(&np_chart, _("Expected Normal"));
1450
1451   chart_initialise(&dnp_chart);
1452   chart_write_title(&dnp_chart, _("Detrended Normal Q-Q Plot of %s"), 
1453                     factorname);
1454   chart_write_xlabel(&dnp_chart, _("Observed Value"));
1455   chart_write_ylabel(&dnp_chart, _("Dev from Normal"));
1456
1457   yfirst = gsl_cdf_ugaussian_Pinv (wv[0].rank / ( m->n + 1));
1458   ylast =  gsl_cdf_ugaussian_Pinv (wv[m->n_data-1].rank / ( m->n + 1));
1459
1460   {
1461     /* Need to make sure that both the scatter plot and the ideal fit into the
1462        plot */
1463     double x_lower = min(m->min, (yfirst - intercept) / slope) ;
1464     double x_upper = max(m->max, (ylast  - intercept) / slope) ;
1465     double slack = (x_upper - x_lower)  * 0.05 ;
1466
1467     chart_write_xscale(&np_chart, x_lower  - slack, x_upper + slack,
1468                        chart_rounded_tick((m->max - m->min) / 5.0));
1469
1470
1471     chart_write_xscale(&dnp_chart, m->min, m->max,
1472                        chart_rounded_tick((m->max - m->min) / 5.0));
1473
1474   }
1475
1476   chart_write_yscale(&np_chart, yfirst, ylast, 
1477                      chart_rounded_tick((ylast - yfirst)/5.0) );
1478
1479   {
1480   /* We have to cache the detrended data, beacause we need to 
1481      find its limits before we can plot it */
1482   double *d_data;
1483   d_data = xmalloc (m->n_data * sizeof(double));
1484   double d_max = -DBL_MAX;
1485   double d_min = DBL_MAX;
1486   for ( i = 0 ; i < m->n_data; ++i ) 
1487     {
1488       const double ns = gsl_cdf_ugaussian_Pinv (wv[i].rank / ( m->n + 1));
1489
1490       chart_datum(&np_chart, 0, wv[i].v.f, ns);
1491
1492       d_data[i] = (wv[i].v.f - m->mean) / m->stddev  - ns;
1493    
1494       if ( d_data[i] < d_min ) d_min = d_data[i];
1495       if ( d_data[i] > d_max ) d_max = d_data[i];
1496     }
1497
1498   chart_write_yscale(&dnp_chart, d_min, d_max, 
1499                      chart_rounded_tick((d_max - d_min) / 5.0));
1500
1501   for ( i = 0 ; i < m->n_data; ++i ) 
1502       chart_datum(&dnp_chart, 0, wv[i].v.f, d_data[i]);
1503
1504   free(d_data);
1505   }
1506
1507   chart_line(&np_chart, slope, intercept, yfirst, ylast , CHART_DIM_Y);
1508   chart_line(&dnp_chart, 0, 0, m->min, m->max , CHART_DIM_X);
1509
1510   chart_finalise(&np_chart);
1511   chart_finalise(&dnp_chart);
1512
1513 }