312547f3f67132380f33d9a204d7049cfb1a670f
[pspp-builds.git] / src / language / stats / examine.q
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2004, 2008 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 <libpspp/message.h>
21 #include <math.h>
22 #include <stdio.h>
23 #include <stdlib.h>
24
25 #include <math/sort.h>
26 #include <math/order-stats.h>
27 #include <math/percentiles.h>
28 #include <math/tukey-hinges.h>
29 #include <math/box-whisker.h>
30 #include <math/trimmed-mean.h>
31 #include <math/extrema.h>
32 #include <math/np.h>
33 #include <data/case.h>
34 #include <data/casegrouper.h>
35 #include <data/casereader.h>
36 #include <data/casewriter.h>
37 #include <data/dictionary.h>
38 #include <data/procedure.h>
39 #include <data/subcase.h>
40 #include <data/value-labels.h>
41 #include <data/variable.h>
42 #include <language/command.h>
43 #include <language/dictionary/split-file.h>
44 #include <language/lexer/lexer.h>
45 #include <libpspp/compiler.h>
46 #include <libpspp/hash.h>
47 #include <libpspp/message.h>
48 #include <libpspp/misc.h>
49 #include <libpspp/str.h>
50 #include <math/moments.h>
51 #include <output/charts/box-whisker.h>
52 #include <output/charts/cartesian.h>
53 #include <output/manager.h>
54 #include <output/table.h>
55
56 #include "minmax.h"
57 #include "xalloc.h"
58
59 #include "gettext.h"
60 #define _(msgid) gettext (msgid)
61 #define N_(msgid) msgid
62
63 /* (headers) */
64 #include <output/chart.h>
65 #include <output/charts/plot-hist.h>
66 #include <output/charts/plot-chart.h>
67 #include <math/histogram.h>
68
69 /* (specification)
70    "EXAMINE" (xmn_):
71    *^variables=custom;
72    +total=custom;
73    +nototal=custom;
74    missing=miss:pairwise/!listwise,
75    rep:report/!noreport,
76    incl:include/!exclude;
77    +compare=cmp:variables/!groups;
78    +percentiles=custom;
79    +id=var;
80    +plot[plt_]=stemleaf,boxplot,npplot,:spreadlevel(*d:n),histogram,all,none;
81    +cinterval=double;
82    +statistics[st_]=descriptives,:extreme(*d:n),all,none.
83 */
84
85 /* (declarations) */
86
87 /* (functions) */
88
89
90 static struct cmd_examine cmd;
91
92 static const struct variable **dependent_vars;
93 static size_t n_dependent_vars;
94
95 /* PERCENTILES */
96
97 static subc_list_double percentile_list;
98 static enum pc_alg percentile_algorithm;
99
100 struct factor_metrics
101 {
102   struct moments1 *moments;
103
104   struct percentile **ptl;
105   size_t n_ptiles;
106
107   struct statistic *tukey_hinges;
108   struct statistic *box_whisker;
109   struct statistic *trimmed_mean;
110   struct statistic *histogram;
111   struct order_stats *np;
112
113   /* Three quartiles indexing into PTL */
114   struct percentile **quartiles;
115
116   /* A reader sorted in ASCENDING order */
117   struct casereader *up_reader;
118
119   /* The minimum value of all the weights */
120   double cmin;
121
122   /* Sum of all weights, including those for missing values */
123   double n;
124
125   double mean;
126
127   double variance;
128
129   double skewness;
130
131   double kurtosis;
132
133   double se_mean;
134
135   struct extrema *minima;
136   struct extrema *maxima;
137 };
138
139 struct factor_result
140 {
141   struct ll ll;
142
143   union value *value[2];
144
145   /* An array of factor metrics, one for each variable */
146   struct factor_metrics *metrics;
147 };
148
149 struct xfactor
150 {
151   /* We need to make a list of this structure */
152   struct ll ll;
153
154   /* The independent variable */
155   const struct variable const* indep_var[2];
156
157   /* A list of results for this factor */
158   struct ll_list result_list ;
159 };
160
161
162 static void
163 factor_destroy (struct xfactor *fctr)
164 {
165   struct ll *ll = ll_head (&fctr->result_list);
166   while (ll != ll_null (&fctr->result_list))
167     {
168       int v;
169       struct factor_result *result =
170         ll_data (ll, struct factor_result, ll);
171
172       for (v = 0; v < n_dependent_vars; ++v)
173         {
174           int i;
175           moments1_destroy (result->metrics[v].moments);
176           extrema_destroy (result->metrics[v].minima);
177           extrema_destroy (result->metrics[v].maxima);
178           statistic_destroy (result->metrics[v].trimmed_mean);
179           statistic_destroy (result->metrics[v].tukey_hinges);
180           statistic_destroy (result->metrics[v].box_whisker);
181           statistic_destroy (result->metrics[v].histogram);
182           for (i = 0 ; i < result->metrics[v].n_ptiles; ++i)
183             statistic_destroy ((struct statistic *) result->metrics[v].ptl[i]);
184           free (result->metrics[v].ptl);
185           free (result->metrics[v].quartiles);
186           casereader_destroy (result->metrics[v].up_reader);
187         }
188
189       free (result->value[0]);
190       free (result->value[1]);
191       free (result->metrics);
192       ll = ll_next (ll);
193       free (result);
194     }
195 }
196
197 static struct xfactor level0_factor;
198 static struct ll_list factor_list = LL_INITIALIZER (factor_list);
199
200 /* Parse the clause specifying the factors */
201 static int examine_parse_independent_vars (struct lexer *lexer,
202                                            const struct dictionary *dict,
203                                            struct cmd_examine *cmd);
204
205 /* Output functions */
206 static void show_summary (const struct variable **dependent_var, int n_dep_var,
207                           const struct xfactor *f);
208
209
210 static void show_descriptives (const struct variable **dependent_var,
211                                int n_dep_var,
212                                const struct xfactor *f);
213
214
215 static void show_percentiles (const struct variable **dependent_var,
216                                int n_dep_var,
217                                const struct xfactor *f);
218
219
220 static void show_extremes (const struct variable **dependent_var,
221                            int n_dep_var,
222                            const struct xfactor *f);
223
224
225
226
227 /* Per Split function */
228 static void run_examine (struct cmd_examine *, struct casereader *,
229                          struct dataset *);
230
231 static void output_examine (void);
232
233
234 void factor_calc (const struct ccase *c, int case_no,
235                   double weight, bool case_missing);
236
237
238 /* Represent a factor as a string, so it can be
239    printed in a human readable fashion */
240 static void factor_to_string (const struct xfactor *fctr,
241                               const struct factor_result *result,
242                               struct string *str);
243
244 /* Represent a factor as a string, so it can be
245    printed in a human readable fashion,
246    but sacrificing some readablility for the sake of brevity */
247 static void
248 factor_to_string_concise (const struct xfactor *fctr,
249                           const struct factor_result *result,
250                           struct string *str
251                           );
252
253
254
255 /* Categories of missing values to exclude. */
256 static enum mv_class exclude_values;
257
258 int
259 cmd_examine (struct lexer *lexer, struct dataset *ds)
260 {
261   struct casegrouper *grouper;
262   struct casereader *group;
263   bool ok;
264
265   subc_list_double_create (&percentile_list);
266   percentile_algorithm = PC_HAVERAGE;
267
268   if ( !parse_examine (lexer, ds, &cmd, NULL) )
269     {
270       subc_list_double_destroy (&percentile_list);
271       return CMD_FAILURE;
272     }
273
274   /* If /MISSING=INCLUDE is set, then user missing values are ignored */
275   exclude_values = cmd.incl == XMN_INCLUDE ? MV_SYSTEM : MV_ANY;
276
277   if ( cmd.st_n == SYSMIS )
278     cmd.st_n = 5;
279
280   if ( ! cmd.sbc_cinterval)
281     cmd.n_cinterval[0] = 95.0;
282
283   /* If descriptives have been requested, make sure the
284      quartiles are calculated */
285   if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
286     {
287       subc_list_double_push (&percentile_list, 25);
288       subc_list_double_push (&percentile_list, 50);
289       subc_list_double_push (&percentile_list, 75);
290     }
291
292   grouper = casegrouper_create_splits (proc_open (ds), dataset_dict (ds));
293
294   while (casegrouper_get_next_group (grouper, &group))
295     {
296       struct casereader *reader =
297         casereader_create_arithmetic_sequence (group, 1, 1);
298
299       run_examine (&cmd, reader, ds);
300     }
301
302   ok = casegrouper_destroy (grouper);
303   ok = proc_commit (ds) && ok;
304
305   if ( dependent_vars )
306     free (dependent_vars);
307
308   subc_list_double_destroy (&percentile_list);
309
310   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
311 };
312
313
314 /* Plot the normal and detrended normal plots for RESULT.
315    Label the plots with LABEL */
316 static void
317 np_plot (struct np *np, const char *label)
318 {
319   double yfirst = 0, ylast = 0;
320
321   double x_lower;
322   double x_upper;
323   double slack;
324
325   /* Normal Plot */
326   struct chart *np_chart;
327
328   /* Detrended Normal Plot */
329   struct chart *dnp_chart;
330
331   /* The slope and intercept of the ideal normal probability line */
332   const double slope = 1.0 / np->stddev;
333   const double intercept = -np->mean / np->stddev;
334
335   if ( np->n < 1.0 )
336     {
337       msg (MW, _("Not creating plot because data set is empty."));
338       return ;
339     }
340
341   np_chart = chart_create ();
342   dnp_chart = chart_create ();
343
344   if ( !np_chart || ! dnp_chart )
345     return ;
346
347   chart_write_title (np_chart, _("Normal Q-Q Plot of %s"), label);
348   chart_write_xlabel (np_chart, _("Observed Value"));
349   chart_write_ylabel (np_chart, _("Expected Normal"));
350
351   chart_write_title (dnp_chart, _("Detrended Normal Q-Q Plot of %s"),
352                      label);
353   chart_write_xlabel (dnp_chart, _("Observed Value"));
354   chart_write_ylabel (dnp_chart, _("Dev from Normal"));
355
356   yfirst = gsl_cdf_ugaussian_Pinv (1 / (np->n + 1));
357   ylast = gsl_cdf_ugaussian_Pinv (np->n / (np->n + 1));
358
359   /* Need to make sure that both the scatter plot and the ideal fit into the
360      plot */
361   x_lower = MIN (np->y_min, (yfirst - intercept) / slope) ;
362   x_upper = MAX (np->y_max, (ylast  - intercept) / slope) ;
363   slack = (x_upper - x_lower)  * 0.05 ;
364
365   chart_write_xscale (np_chart, x_lower - slack, x_upper + slack, 5);
366   chart_write_xscale (dnp_chart, np->y_min, np->y_max, 5);
367
368   chart_write_yscale (np_chart, yfirst, ylast, 5);
369   chart_write_yscale (dnp_chart, np->dns_min, np->dns_max, 5);
370
371   {
372     struct ccase c;
373     struct casereader *reader = casewriter_make_reader (np->writer);
374     while (casereader_read (reader, &c))
375       {
376         chart_datum (np_chart, 0, case_data_idx (&c, NP_IDX_Y)->f, case_data_idx (&c, NP_IDX_NS)->f);
377         chart_datum (dnp_chart, 0, case_data_idx (&c, NP_IDX_Y)->f, case_data_idx (&c, NP_IDX_DNS)->f);
378
379         case_destroy (&c);
380       }
381     casereader_destroy (reader);
382   }
383
384   chart_line (dnp_chart, 0, 0, np->y_min, np->y_max , CHART_DIM_X);
385   chart_line (np_chart, slope, intercept, yfirst, ylast , CHART_DIM_Y);
386
387   chart_submit (np_chart);
388   chart_submit (dnp_chart);
389 }
390
391
392 static void
393 show_npplot (const struct variable **dependent_var,
394              int n_dep_var,
395              const struct xfactor *fctr)
396 {
397   int v;
398
399   for (v = 0; v < n_dep_var; ++v)
400     {
401       struct ll *ll;
402       for (ll = ll_head (&fctr->result_list);
403            ll != ll_null (&fctr->result_list);
404            ll = ll_next (ll))
405         {
406           struct string str;
407           const struct factor_result *result =
408             ll_data (ll, struct factor_result, ll);
409
410           ds_init_empty (&str);
411           ds_put_format (&str, "%s ", var_get_name (dependent_var[v]));
412
413           factor_to_string (fctr, result, &str);
414
415           np_plot ((struct np*) result->metrics[v].np, ds_cstr(&str));
416
417           statistic_destroy ((struct statistic *)result->metrics[v].np);
418
419           ds_destroy (&str);
420         }
421     }
422 }
423
424
425 static void
426 show_histogram (const struct variable **dependent_var,
427                 int n_dep_var,
428                 const struct xfactor *fctr)
429 {
430   int v;
431
432   for (v = 0; v < n_dep_var; ++v)
433     {
434       struct ll *ll;
435       for (ll = ll_head (&fctr->result_list);
436            ll != ll_null (&fctr->result_list);
437            ll = ll_next (ll))
438         {
439           struct string str;
440           const struct factor_result *result =
441             ll_data (ll, struct factor_result, ll);
442
443           ds_init_empty (&str);
444           ds_put_format (&str, "%s ", var_get_name (dependent_var[v]));
445
446           factor_to_string (fctr, result, &str);
447
448           histogram_plot ((struct histogram *) result->metrics[v].histogram,
449                           ds_cstr (&str),
450                           (struct moments1 *) result->metrics[v].moments);
451
452           ds_destroy (&str);
453         }
454     }
455 }
456
457
458
459 static void
460 show_boxplot_groups (const struct variable **dependent_var,
461                      int n_dep_var,
462                      const struct xfactor *fctr)
463 {
464   int v;
465
466   for (v = 0; v < n_dep_var; ++v)
467     {
468       struct ll *ll;
469       int f = 0;
470       struct chart *ch = chart_create ();
471       double y_min = DBL_MAX;
472       double y_max = -DBL_MAX;
473
474       for (ll = ll_head (&fctr->result_list);
475            ll != ll_null (&fctr->result_list);
476            ll = ll_next (ll))
477         {
478           const struct extremum  *max, *min;
479           const struct factor_result *result =
480             ll_data (ll, struct factor_result, ll);
481
482           const struct ll_list *max_list =
483             extrema_list (result->metrics[v].maxima);
484
485           const struct ll_list *min_list =
486             extrema_list (result->metrics[v].minima);
487
488           if ( ll_is_empty (max_list))
489             {
490               msg (MW, _("Not creating plot because data set is empty."));
491               continue;
492             }
493
494           max = (const struct extremum *)
495             ll_data (ll_head(max_list), struct extremum, ll);
496
497           min = (const struct extremum *)
498             ll_data (ll_head (min_list), struct extremum, ll);
499
500           y_max = MAX (y_max, max->value);
501           y_min = MIN (y_min, min->value);
502         }
503
504       boxplot_draw_yscale (ch, y_max, y_min);
505
506       if ( fctr->indep_var[0])
507         chart_write_title (ch, _("Boxplot of %s vs. %s"),
508                            var_to_string (dependent_var[v]),
509                            var_to_string (fctr->indep_var[0]) );
510       else
511         chart_write_title (ch, _("Boxplot of %s"),
512                            var_to_string (dependent_var[v]));
513
514       for (ll = ll_head (&fctr->result_list);
515            ll != ll_null (&fctr->result_list);
516            ll = ll_next (ll))
517         {
518           const struct factor_result *result =
519             ll_data (ll, struct factor_result, ll);
520
521           struct string str;
522           const double box_width = (ch->data_right - ch->data_left)
523             / (ll_count (&fctr->result_list) * 2.0 ) ;
524
525           const double box_centre = (f++ * 2 + 1) * box_width + ch->data_left;
526
527           ds_init_empty (&str);
528           factor_to_string_concise (fctr, result, &str);
529
530           boxplot_draw_boxplot (ch,
531                                 box_centre, box_width,
532                                 (const struct box_whisker *)
533                                  result->metrics[v].box_whisker,
534                                 ds_cstr (&str));
535
536           ds_destroy (&str);
537         }
538
539       chart_submit (ch);
540     }
541 }
542
543
544
545 static void
546 show_boxplot_variables (const struct variable **dependent_var,
547                         int n_dep_var,
548                         const struct xfactor *fctr
549                         )
550
551 {
552   int v;
553   struct ll *ll;
554   const struct ll_list *result_list = &fctr->result_list;
555
556   for (ll = ll_head (result_list);
557        ll != ll_null (result_list);
558        ll = ll_next (ll))
559
560     {
561       struct string title;
562       struct chart *ch = chart_create ();
563       double y_min = DBL_MAX;
564       double y_max = -DBL_MAX;
565
566       const struct factor_result *result =
567         ll_data (ll, struct factor_result, ll);
568
569       const double box_width = (ch->data_right - ch->data_left)
570         / (n_dep_var * 2.0 ) ;
571
572       for (v = 0; v < n_dep_var; ++v)
573         {
574           const struct ll *max_ll =
575             ll_head (extrema_list (result->metrics[v].maxima));
576           const struct ll *min_ll =
577             ll_head (extrema_list (result->metrics[v].minima));
578
579           const struct extremum  *max =
580             (const struct extremum *) ll_data (max_ll, struct extremum, ll);
581
582           const struct extremum  *min =
583             (const struct extremum *) ll_data (min_ll, struct extremum, ll);
584
585           y_max = MAX (y_max, max->value);
586           y_min = MIN (y_min, min->value);
587         }
588
589
590       boxplot_draw_yscale (ch, y_max, y_min);
591
592       ds_init_empty (&title);
593       factor_to_string (fctr, result, &title);
594
595 #if 0
596       ds_put_format (&title, "%s = ", var_get_name (fctr->indep_var[0]));
597       var_append_value_name (fctr->indep_var[0], result->value[0], &title);
598 #endif
599
600       chart_write_title (ch, ds_cstr (&title));
601       ds_destroy (&title);
602
603       for (v = 0; v < n_dep_var; ++v)
604         {
605           struct string str;
606           const double box_centre = (v * 2 + 1) * box_width + ch->data_left;
607
608           ds_init_empty (&str);
609           ds_init_cstr (&str, var_get_name (dependent_var[v]));
610
611           boxplot_draw_boxplot (ch,
612                                 box_centre, box_width,
613                                 (const struct box_whisker *) result->metrics[v].box_whisker,
614                                 ds_cstr (&str));
615
616           ds_destroy (&str);
617         }
618
619       chart_submit (ch);
620     }
621 }
622
623
624 /* Show all the appropriate tables */
625 static void
626 output_examine (void)
627 {
628   struct ll *ll;
629
630   show_summary (dependent_vars, n_dependent_vars, &level0_factor);
631
632   if ( cmd.a_statistics[XMN_ST_EXTREME] )
633     show_extremes (dependent_vars, n_dependent_vars, &level0_factor);
634
635   if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
636     show_descriptives (dependent_vars, n_dependent_vars, &level0_factor);
637
638   if ( cmd.sbc_percentiles)
639     show_percentiles (dependent_vars, n_dependent_vars, &level0_factor);
640
641   if ( cmd.sbc_plot)
642     {
643       if (cmd.a_plot[XMN_PLT_BOXPLOT])
644         show_boxplot_groups (dependent_vars, n_dependent_vars, &level0_factor);
645
646       if (cmd.a_plot[XMN_PLT_HISTOGRAM])
647         show_histogram (dependent_vars, n_dependent_vars, &level0_factor);
648
649       if (cmd.a_plot[XMN_PLT_NPPLOT])
650         show_npplot (dependent_vars, n_dependent_vars, &level0_factor);
651     }
652
653   for (ll = ll_head (&factor_list);
654        ll != ll_null (&factor_list); ll = ll_next (ll))
655     {
656       struct xfactor *factor = ll_data (ll, struct xfactor, ll);
657       show_summary (dependent_vars, n_dependent_vars, factor);
658
659       if ( cmd.a_statistics[XMN_ST_EXTREME] )
660         show_extremes (dependent_vars, n_dependent_vars, factor);
661
662       if ( cmd.a_statistics[XMN_ST_DESCRIPTIVES] )
663         show_descriptives (dependent_vars, n_dependent_vars, factor);
664
665       if ( cmd.sbc_percentiles)
666         show_percentiles (dependent_vars, n_dependent_vars, factor);
667
668       if (cmd.a_plot[XMN_PLT_BOXPLOT] &&
669           cmd.cmp == XMN_GROUPS)
670         show_boxplot_groups (dependent_vars, n_dependent_vars, factor);
671
672
673       if (cmd.a_plot[XMN_PLT_BOXPLOT] &&
674           cmd.cmp == XMN_VARIABLES)
675         show_boxplot_variables (dependent_vars, n_dependent_vars,
676                                 factor);
677
678       if (cmd.a_plot[XMN_PLT_HISTOGRAM])
679         show_histogram (dependent_vars, n_dependent_vars, factor);
680
681       if (cmd.a_plot[XMN_PLT_NPPLOT])
682         show_npplot (dependent_vars, n_dependent_vars, factor);
683     }
684 }
685
686 /* Parse the PERCENTILES subcommand */
687 static int
688 xmn_custom_percentiles (struct lexer *lexer, struct dataset *ds UNUSED,
689                         struct cmd_examine *p UNUSED, void *aux UNUSED)
690 {
691   lex_match (lexer, '=');
692
693   lex_match (lexer, '(');
694
695   while ( lex_is_number (lexer) )
696     {
697       subc_list_double_push (&percentile_list, lex_number (lexer));
698
699       lex_get (lexer);
700
701       lex_match (lexer, ',') ;
702     }
703   lex_match (lexer, ')');
704
705   lex_match (lexer, '=');
706
707   if ( lex_match_id (lexer, "HAVERAGE"))
708     percentile_algorithm = PC_HAVERAGE;
709
710   else if ( lex_match_id (lexer, "WAVERAGE"))
711     percentile_algorithm = PC_WAVERAGE;
712
713   else if ( lex_match_id (lexer, "ROUND"))
714     percentile_algorithm = PC_ROUND;
715
716   else if ( lex_match_id (lexer, "EMPIRICAL"))
717     percentile_algorithm = PC_EMPIRICAL;
718
719   else if ( lex_match_id (lexer, "AEMPIRICAL"))
720     percentile_algorithm = PC_AEMPIRICAL;
721
722   else if ( lex_match_id (lexer, "NONE"))
723     percentile_algorithm = PC_NONE;
724
725
726   if ( 0 == subc_list_double_count (&percentile_list))
727     {
728       subc_list_double_push (&percentile_list, 5);
729       subc_list_double_push (&percentile_list, 10);
730       subc_list_double_push (&percentile_list, 25);
731       subc_list_double_push (&percentile_list, 50);
732       subc_list_double_push (&percentile_list, 75);
733       subc_list_double_push (&percentile_list, 90);
734       subc_list_double_push (&percentile_list, 95);
735     }
736
737   return 1;
738 }
739
740 /* TOTAL and NOTOTAL are simple, mutually exclusive flags */
741 static int
742 xmn_custom_total (struct lexer *lexer UNUSED, struct dataset *ds UNUSED,
743                   struct cmd_examine *p, void *aux UNUSED)
744 {
745   if ( p->sbc_nototal )
746     {
747       msg (SE, _("%s and %s are mutually exclusive"),"TOTAL","NOTOTAL");
748       return 0;
749     }
750
751   return 1;
752 }
753
754 static int
755 xmn_custom_nototal (struct lexer *lexer UNUSED, struct dataset *ds UNUSED,
756                     struct cmd_examine *p, void *aux UNUSED)
757 {
758   if ( p->sbc_total )
759     {
760       msg (SE, _("%s and %s are mutually exclusive"), "TOTAL", "NOTOTAL");
761       return 0;
762     }
763
764   return 1;
765 }
766
767
768
769 /* Parser for the variables sub command
770    Returns 1 on success */
771 static int
772 xmn_custom_variables (struct lexer *lexer, struct dataset *ds,
773                       struct cmd_examine *cmd,
774                       void *aux UNUSED)
775 {
776   const struct dictionary *dict = dataset_dict (ds);
777   lex_match (lexer, '=');
778
779   if ( (lex_token (lexer) != T_ID || dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
780        && lex_token (lexer) != T_ALL)
781     {
782       return 2;
783     }
784
785   if (!parse_variables_const (lexer, dict, &dependent_vars, &n_dependent_vars,
786                               PV_NO_DUPLICATE | PV_NUMERIC | PV_NO_SCRATCH) )
787     {
788       free (dependent_vars);
789       return 0;
790     }
791
792   assert (n_dependent_vars);
793
794
795   if ( lex_match (lexer, T_BY))
796     {
797       int success ;
798       success =  examine_parse_independent_vars (lexer, dict, cmd);
799       if ( success != 1 )
800         {
801           free (dependent_vars);
802         }
803       return success;
804     }
805
806   return 1;
807 }
808
809
810
811 /* Parse the clause specifying the factors */
812 static int
813 examine_parse_independent_vars (struct lexer *lexer,
814                                 const struct dictionary *dict,
815                                 struct cmd_examine *cmd)
816 {
817   int success;
818   struct xfactor *sf = xmalloc (sizeof *sf);
819
820   ll_init (&sf->result_list);
821
822   if ( (lex_token (lexer) != T_ID ||
823         dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
824        && lex_token (lexer) != T_ALL)
825     {
826       free ( sf ) ;
827       return 2;
828     }
829
830   sf->indep_var[0] = parse_variable (lexer, dict);
831   sf->indep_var[1] = NULL;
832
833   if ( lex_token (lexer) == T_BY )
834     {
835       lex_match (lexer, T_BY);
836
837       if ( (lex_token (lexer) != T_ID ||
838             dict_lookup_var (dict, lex_tokid (lexer)) == NULL)
839            && lex_token (lexer) != T_ALL)
840         {
841           free (sf);
842           return 2;
843         }
844
845       sf->indep_var[1] = parse_variable (lexer, dict);
846
847       ll_push_tail (&factor_list, &sf->ll);
848     }
849   else
850     ll_push_tail (&factor_list, &sf->ll);
851
852   lex_match (lexer, ',');
853
854   if ( lex_token (lexer) == '.' || lex_token (lexer) == '/' )
855     return 1;
856
857   success =  examine_parse_independent_vars (lexer, dict, cmd);
858
859   if ( success != 1 )
860     free ( sf ) ;
861
862   return success;
863 }
864
865 static void
866 examine_group (struct cmd_examine *cmd, struct casereader *reader, int level,
867                const struct dictionary *dict, struct xfactor *factor)
868 {
869   struct ccase c;
870   const struct variable *wv = dict_get_weight (dict);
871   int v;
872   int n_extrema = 1;
873   struct factor_result *result = xzalloc (sizeof (*result));
874
875   result->metrics = xcalloc (n_dependent_vars, sizeof (*result->metrics));
876
877   if ( cmd->a_statistics[XMN_ST_EXTREME] )
878     n_extrema = cmd->st_n;
879
880
881   if (casereader_peek (reader, 0, &c))
882     {
883       if ( level > 0)
884         {
885           result->value[0] =
886             value_dup (case_data (&c, factor->indep_var[0]),
887                        var_get_width (factor->indep_var[0]));
888
889           if ( level > 1)
890             result->value[1] =
891               value_dup (case_data (&c, factor->indep_var[1]),
892                          var_get_width (factor->indep_var[1]));
893         }
894       case_destroy (&c);
895     }
896
897   for (v = 0; v < n_dependent_vars; ++v)
898     {
899       struct casewriter *writer;
900       struct casereader *input = casereader_clone (reader);
901
902       result->metrics[v].moments = moments1_create (MOMENT_KURTOSIS);
903       result->metrics[v].minima = extrema_create (n_extrema, EXTREME_MINIMA);
904       result->metrics[v].maxima = extrema_create (n_extrema, EXTREME_MAXIMA);
905       result->metrics[v].cmin = DBL_MAX;
906
907       if (cmd->a_statistics[XMN_ST_DESCRIPTIVES] ||
908           cmd->a_plot[XMN_PLT_BOXPLOT] ||
909           cmd->a_plot[XMN_PLT_NPPLOT] ||
910           cmd->sbc_percentiles)
911         {
912           /* In this case, we need to sort the data, so we create a sorting
913              casewriter */
914           struct subcase up_ordering;
915           subcase_init_var (&up_ordering, dependent_vars[v], SC_ASCEND);
916           writer = sort_create_writer (&up_ordering,
917                                        casereader_get_value_cnt (reader));
918           subcase_destroy (&up_ordering);
919         }
920       else
921         {
922           /* but in this case, sorting is unnecessary, so an ordinary
923              casewriter is sufficient */
924           writer =
925             autopaging_writer_create (casereader_get_value_cnt (reader));
926         }
927
928
929       /* Sort or just iterate, whilst calculating moments etc */
930       while (casereader_read (input, &c))
931         {
932           const casenumber loc =
933               case_data_idx (&c, casereader_get_value_cnt (reader) - 1)->f;
934
935           const double weight = wv ? case_data (&c, wv)->f : 1.0;
936
937           if (weight != SYSMIS)
938             minimize (&result->metrics[v].cmin, weight);
939
940           moments1_add (result->metrics[v].moments,
941                         case_data (&c, dependent_vars[v])->f,
942                         weight);
943
944           result->metrics[v].n += weight;
945
946           extrema_add (result->metrics[v].maxima,
947                        case_data (&c, dependent_vars[v])->f,
948                        weight,
949                        loc);
950
951           extrema_add (result->metrics[v].minima,
952                        case_data (&c, dependent_vars[v])->f,
953                        weight,
954                        loc);
955
956           casewriter_write (writer, &c);
957         }
958       casereader_destroy (input);
959       result->metrics[v].up_reader = casewriter_make_reader (writer);
960     }
961
962   /* If percentiles or descriptives have been requested, then a
963      second pass through the data (which has now been sorted)
964      is necessary */
965   if ( cmd->a_statistics[XMN_ST_DESCRIPTIVES] ||
966        cmd->a_plot[XMN_PLT_BOXPLOT] ||
967        cmd->a_plot[XMN_PLT_NPPLOT] ||
968        cmd->sbc_percentiles)
969     {
970       for (v = 0; v < n_dependent_vars; ++v)
971         {
972           int i;
973           int n_os;
974           struct order_stats **os ;
975           struct factor_metrics *metric = &result->metrics[v];
976
977           metric->n_ptiles = percentile_list.n_data;
978
979           metric->ptl = xcalloc (metric->n_ptiles,
980                                  sizeof (struct percentile *));
981
982           metric->quartiles = xcalloc (3, sizeof (*metric->quartiles));
983
984           for (i = 0 ; i < metric->n_ptiles; ++i)
985             {
986               metric->ptl[i] = (struct percentile *)
987                 percentile_create (percentile_list.data[i] / 100.0, metric->n);
988
989               if ( percentile_list.data[i] == 25)
990                 metric->quartiles[0] = metric->ptl[i];
991               else if ( percentile_list.data[i] == 50)
992                 metric->quartiles[1] = metric->ptl[i];
993               else if ( percentile_list.data[i] == 75)
994                 metric->quartiles[2] = metric->ptl[i];
995             }
996
997           metric->tukey_hinges = tukey_hinges_create (metric->n, metric->cmin);
998           metric->trimmed_mean = trimmed_mean_create (metric->n, 0.05);
999
1000           n_os = metric->n_ptiles + 2;
1001
1002          if ( cmd->a_plot[XMN_PLT_NPPLOT] )
1003             {
1004               metric->np = np_create (metric->moments);
1005               n_os ++;
1006             }
1007
1008           os = xcalloc (sizeof (struct order_stats *), n_os);
1009
1010           for (i = 0 ; i < metric->n_ptiles ; ++i )
1011             {
1012               os[i] = (struct order_stats *) metric->ptl[i];
1013             }
1014
1015           os[i] = (struct order_stats *) metric->tukey_hinges;
1016           os[i+1] = (struct order_stats *) metric->trimmed_mean;
1017
1018           if (cmd->a_plot[XMN_PLT_NPPLOT])
1019             os[i+2] = metric->np;
1020
1021           order_stats_accumulate (os, n_os,
1022                                   casereader_clone (metric->up_reader),
1023                                   wv, dependent_vars[v], MV_ANY);
1024           free (os);
1025         }
1026     }
1027
1028   /* FIXME: Do this in the above loop */
1029   if ( cmd->a_plot[XMN_PLT_HISTOGRAM] )
1030     {
1031       struct ccase c;
1032       struct casereader *input = casereader_clone (reader);
1033
1034       for (v = 0; v < n_dependent_vars; ++v)
1035         {
1036           const struct extremum  *max, *min;
1037           struct factor_metrics *metric = &result->metrics[v];
1038
1039           const struct ll_list *max_list =
1040             extrema_list (result->metrics[v].maxima);
1041
1042           const struct ll_list *min_list =
1043             extrema_list (result->metrics[v].minima);
1044
1045           if ( ll_is_empty (max_list))
1046             {
1047               msg (MW, _("Not creating plot because data set is empty."));
1048               continue;
1049             }
1050
1051           assert (! ll_is_empty (min_list));
1052
1053           max = (const struct extremum *)
1054             ll_data (ll_head(max_list), struct extremum, ll);
1055
1056           min = (const struct extremum *)
1057             ll_data (ll_head (min_list), struct extremum, ll);
1058
1059           metric->histogram = histogram_create (10, min->value, max->value);
1060         }
1061
1062       while (casereader_read (input, &c))
1063         {
1064           const double weight = wv ? case_data (&c, wv)->f : 1.0;
1065
1066           for (v = 0; v < n_dependent_vars; ++v)
1067             {
1068               struct factor_metrics *metric = &result->metrics[v];
1069               if ( metric->histogram)
1070                 histogram_add ((struct histogram *) metric->histogram,
1071                                case_data (&c, dependent_vars[v])->f, weight);
1072             }
1073           case_destroy (&c);
1074         }
1075       casereader_destroy (input);
1076     }
1077
1078   /* In this case, a third iteration is required */
1079   if (cmd->a_plot[XMN_PLT_BOXPLOT])
1080     {
1081       for (v = 0; v < n_dependent_vars; ++v)
1082         {
1083           struct factor_metrics *metric = &result->metrics[v];
1084
1085           metric->box_whisker =
1086             box_whisker_create ((struct tukey_hinges *) metric->tukey_hinges,
1087                                 cmd->v_id,
1088                                 casereader_get_value_cnt (metric->up_reader)
1089                                 - 1);
1090
1091           order_stats_accumulate ((struct order_stats **) &metric->box_whisker,
1092                                   1,
1093                                   casereader_clone (metric->up_reader),
1094                                   wv, dependent_vars[v], MV_ANY);
1095         }
1096     }
1097
1098   ll_push_tail (&factor->result_list, &result->ll);
1099   casereader_destroy (reader);
1100 }
1101
1102
1103 static void
1104 run_examine (struct cmd_examine *cmd, struct casereader *input,
1105              struct dataset *ds)
1106 {
1107   struct ll *ll;
1108   const struct dictionary *dict = dataset_dict (ds);
1109   struct ccase c;
1110   struct casereader *level0 = casereader_clone (input);
1111
1112   if (!casereader_peek (input, 0, &c))
1113     {
1114       casereader_destroy (input);
1115       return;
1116     }
1117
1118   output_split_file_values (ds, &c);
1119   case_destroy (&c);
1120
1121   ll_init (&level0_factor.result_list);
1122
1123   examine_group (cmd, level0, 0, dict, &level0_factor);
1124
1125   for (ll = ll_head (&factor_list);
1126        ll != ll_null (&factor_list);
1127        ll = ll_next (ll))
1128     {
1129       struct xfactor *factor = ll_data (ll, struct xfactor, ll);
1130
1131       struct casereader *group = NULL;
1132       struct casereader *level1;
1133       struct casegrouper *grouper1 = NULL;
1134
1135       level1 = casereader_clone (input);
1136       level1 = sort_execute_1var (level1, factor->indep_var[0]);
1137       grouper1 = casegrouper_create_vars (level1, &factor->indep_var[0], 1);
1138
1139       while (casegrouper_get_next_group (grouper1, &group))
1140         {
1141           struct casereader *group_copy = casereader_clone (group);
1142
1143           if ( !factor->indep_var[1])
1144             examine_group (cmd, group_copy, 1, dict, factor);
1145           else
1146             {
1147               int n_groups = 0;
1148               struct casereader *group2 = NULL;
1149               struct casegrouper *grouper2 = NULL;
1150
1151               group_copy = sort_execute_1var (group_copy,
1152                                               factor->indep_var[1]);
1153
1154               grouper2 = casegrouper_create_vars (group_copy,
1155                                                   &factor->indep_var[1], 1);
1156
1157               while (casegrouper_get_next_group (grouper2, &group2))
1158                 {
1159                   examine_group (cmd, group2, 2, dict, factor);
1160                   n_groups++;
1161                 }
1162               casegrouper_destroy (grouper2);
1163             }
1164
1165           casereader_destroy (group);
1166         }
1167       casegrouper_destroy (grouper1);
1168     }
1169
1170   casereader_destroy (input);
1171
1172   output_examine ();
1173
1174   factor_destroy (&level0_factor);
1175
1176   {
1177     struct ll *ll;
1178     for (ll = ll_head (&factor_list);
1179          ll != ll_null (&factor_list);
1180          ll = ll_next (ll))
1181       {
1182         struct xfactor *f = ll_data (ll, struct xfactor, ll);
1183         factor_destroy (f);
1184       }
1185   }
1186
1187 }
1188
1189
1190 static void
1191 show_summary (const struct variable **dependent_var, int n_dep_var,
1192               const struct xfactor *fctr)
1193 {
1194   static const char *subtitle[]=
1195     {
1196       N_("Valid"),
1197       N_("Missing"),
1198       N_("Total")
1199     };
1200
1201   int v, j;
1202   int heading_columns = 1;
1203   int n_cols;
1204   const int heading_rows = 3;
1205   struct tab_table *tbl;
1206
1207   int n_rows ;
1208   n_rows = n_dep_var;
1209
1210   assert (fctr);
1211
1212   if ( fctr->indep_var[0] )
1213     {
1214       heading_columns = 2;
1215
1216       if ( fctr->indep_var[1] )
1217         {
1218           heading_columns = 3;
1219         }
1220     }
1221
1222   n_rows *= ll_count (&fctr->result_list);
1223   n_rows += heading_rows;
1224
1225   n_cols = heading_columns + 6;
1226
1227   tbl = tab_create (n_cols, n_rows, 0);
1228   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1229
1230   tab_dim (tbl, tab_natural_dimensions);
1231
1232   /* Outline the box */
1233   tab_box (tbl,
1234            TAL_2, TAL_2,
1235            -1, -1,
1236            0, 0,
1237            n_cols - 1, n_rows - 1);
1238
1239   /* Vertical lines for the data only */
1240   tab_box (tbl,
1241            -1, -1,
1242            -1, TAL_1,
1243            heading_columns, 0,
1244            n_cols - 1, n_rows - 1);
1245
1246
1247   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1248   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, 1 );
1249   tab_hline (tbl, TAL_1, heading_columns, n_cols - 1, heading_rows -1 );
1250
1251   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
1252
1253
1254   tab_title (tbl, _("Case Processing Summary"));
1255
1256   tab_joint_text (tbl, heading_columns, 0,
1257                   n_cols -1, 0,
1258                   TAB_CENTER | TAT_TITLE,
1259                   _("Cases"));
1260
1261   /* Remove lines ... */
1262   tab_box (tbl,
1263            -1, -1,
1264            TAL_0, TAL_0,
1265            heading_columns, 0,
1266            n_cols - 1, 0);
1267
1268   for (j = 0 ; j < 3 ; ++j)
1269     {
1270       tab_text (tbl, heading_columns + j * 2 , 2, TAB_CENTER | TAT_TITLE,
1271                 _("N"));
1272
1273       tab_text (tbl, heading_columns + j * 2 + 1, 2, TAB_CENTER | TAT_TITLE,
1274                 _("Percent"));
1275
1276       tab_joint_text (tbl, heading_columns + j * 2 , 1,
1277                       heading_columns + j * 2 + 1, 1,
1278                       TAB_CENTER | TAT_TITLE,
1279                       subtitle[j]);
1280
1281       tab_box (tbl, -1, -1,
1282                TAL_0, TAL_0,
1283                heading_columns + j * 2, 1,
1284                heading_columns + j * 2 + 1, 1);
1285     }
1286
1287
1288   /* Titles for the independent variables */
1289   if ( fctr->indep_var[0] )
1290     {
1291       tab_text (tbl, 1, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1292                 var_to_string (fctr->indep_var[0]));
1293
1294       if ( fctr->indep_var[1] )
1295         {
1296           tab_text (tbl, 2, heading_rows - 1, TAB_CENTER | TAT_TITLE,
1297                     var_to_string (fctr->indep_var[1]));
1298         }
1299     }
1300
1301   for (v = 0 ; v < n_dep_var ; ++v)
1302     {
1303       int j = 0;
1304       struct ll *ll;
1305       union value *last_value = NULL;
1306
1307       if ( v > 0 )
1308         tab_hline (tbl, TAL_1, 0, n_cols -1 ,
1309                    v * ll_count (&fctr->result_list)
1310                    + heading_rows);
1311
1312       tab_text (tbl,
1313                 0,
1314                 v * ll_count (&fctr->result_list) + heading_rows,
1315                 TAB_LEFT | TAT_TITLE,
1316                 var_to_string (dependent_var[v])
1317                 );
1318
1319
1320       for (ll = ll_head (&fctr->result_list);
1321            ll != ll_null (&fctr->result_list); ll = ll_next (ll))
1322         {
1323           double n;
1324           const struct factor_result *result =
1325             ll_data (ll, struct factor_result, ll);
1326
1327           if ( fctr->indep_var[0] )
1328             {
1329
1330               if ( last_value == NULL ||
1331                    compare_values_short (last_value, result->value[0],
1332                                          fctr->indep_var[0]))
1333                 {
1334                   struct string str;
1335
1336                   last_value = result->value[0];
1337                   ds_init_empty (&str);
1338
1339                   var_append_value_name (fctr->indep_var[0], result->value[0],
1340                                          &str);
1341
1342                   tab_text (tbl, 1,
1343                             heading_rows + j +
1344                             v * ll_count (&fctr->result_list),
1345                             TAB_LEFT | TAT_TITLE,
1346                             ds_cstr (&str));
1347
1348                   ds_destroy (&str);
1349
1350                   if ( fctr->indep_var[1] && j > 0)
1351                     tab_hline (tbl, TAL_1, 1, n_cols - 1,
1352                                heading_rows + j +
1353                                v * ll_count (&fctr->result_list));
1354                 }
1355
1356               if ( fctr->indep_var[1])
1357                 {
1358                   struct string str;
1359
1360                   ds_init_empty (&str);
1361
1362                   var_append_value_name (fctr->indep_var[1],
1363                                          result->value[1], &str);
1364
1365                   tab_text (tbl, 2,
1366                             heading_rows + j +
1367                             v * ll_count (&fctr->result_list),
1368                             TAB_LEFT | TAT_TITLE,
1369                             ds_cstr (&str));
1370
1371                   ds_destroy (&str);
1372                 }
1373             }
1374
1375
1376           moments1_calculate (result->metrics[v].moments,
1377                               &n, &result->metrics[v].mean,
1378                               &result->metrics[v].variance,
1379                               &result->metrics[v].skewness,
1380                               &result->metrics[v].kurtosis);
1381
1382           result->metrics[v].se_mean = sqrt (result->metrics[v].variance / n) ;
1383
1384           /* Total Valid */
1385           tab_float (tbl, heading_columns,
1386                      heading_rows + j + v * ll_count (&fctr->result_list),
1387                      TAB_LEFT,
1388                      n, 8, 0);
1389
1390           tab_text (tbl, heading_columns + 1,
1391                     heading_rows + j + v * ll_count (&fctr->result_list),
1392                     TAB_RIGHT | TAT_PRINTF,
1393                     "%g%%", n * 100.0 / result->metrics[v].n);
1394
1395           /* Total Missing */
1396           tab_float (tbl, heading_columns + 2,
1397                      heading_rows + j + v * ll_count (&fctr->result_list),
1398                      TAB_LEFT,
1399                      result->metrics[v].n - n,
1400                      8, 0);
1401
1402           tab_text (tbl, heading_columns + 3,
1403                     heading_rows + j + v * ll_count (&fctr->result_list),
1404                     TAB_RIGHT | TAT_PRINTF,
1405                     "%g%%",
1406                     (result->metrics[v].n - n) * 100.0 / result->metrics[v].n
1407                     );
1408
1409           /* Total Valid + Missing */
1410           tab_float (tbl, heading_columns + 4,
1411                      heading_rows + j + v * ll_count (&fctr->result_list),
1412                      TAB_LEFT,
1413                      result->metrics[v].n,
1414                      8, 0);
1415
1416           tab_text (tbl, heading_columns + 5,
1417                     heading_rows + j + v * ll_count (&fctr->result_list),
1418                     TAB_RIGHT | TAT_PRINTF,
1419                     "%g%%",
1420                     (result->metrics[v].n) * 100.0 / result->metrics[v].n
1421                     );
1422
1423           ++j;
1424         }
1425     }
1426
1427
1428   tab_submit (tbl);
1429 }
1430
1431 #define DESCRIPTIVE_ROWS 13
1432
1433 static void
1434 show_descriptives (const struct variable **dependent_var,
1435                    int n_dep_var,
1436                    const struct xfactor *fctr)
1437 {
1438   int v;
1439   int heading_columns = 3;
1440   int n_cols;
1441   const int heading_rows = 1;
1442   struct tab_table *tbl;
1443
1444   int n_rows ;
1445   n_rows = n_dep_var;
1446
1447   assert (fctr);
1448
1449   if ( fctr->indep_var[0] )
1450     {
1451       heading_columns = 4;
1452
1453       if ( fctr->indep_var[1] )
1454         {
1455           heading_columns = 5;
1456         }
1457     }
1458
1459   n_rows *= ll_count (&fctr->result_list) * DESCRIPTIVE_ROWS;
1460   n_rows += heading_rows;
1461
1462   n_cols = heading_columns + 2;
1463
1464   tbl = tab_create (n_cols, n_rows, 0);
1465   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1466
1467   tab_dim (tbl, tab_natural_dimensions);
1468
1469   /* Outline the box */
1470   tab_box (tbl,
1471            TAL_2, TAL_2,
1472            -1, -1,
1473            0, 0,
1474            n_cols - 1, n_rows - 1);
1475
1476
1477   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1478   tab_hline (tbl, TAL_2, 1, n_cols - 1, heading_rows );
1479
1480   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1481
1482
1483   if ( fctr->indep_var[0])
1484     tab_text (tbl, 1, 0, TAT_TITLE, var_to_string (fctr->indep_var[0]));
1485
1486   if ( fctr->indep_var[1])
1487     tab_text (tbl, 2, 0, TAT_TITLE, var_to_string (fctr->indep_var[1]));
1488
1489   for (v = 0 ; v < n_dep_var ; ++v )
1490     {
1491       struct ll *ll;
1492       int i = 0;
1493
1494       const int row_var_start =
1495         v * DESCRIPTIVE_ROWS * ll_count(&fctr->result_list);
1496
1497       tab_text (tbl,
1498                 0,
1499                 heading_rows + row_var_start,
1500                 TAB_LEFT | TAT_TITLE,
1501                 var_to_string (dependent_var[v])
1502                 );
1503
1504       for (ll = ll_head (&fctr->result_list);
1505            ll != ll_null (&fctr->result_list); i++, ll = ll_next (ll))
1506         {
1507           const struct factor_result *result =
1508             ll_data (ll, struct factor_result, ll);
1509
1510           const double t =
1511             gsl_cdf_tdist_Qinv ((1 - cmd.n_cinterval[0] / 100.0) / 2.0,
1512                                       result->metrics[v].n - 1);
1513
1514           if ( i > 0 || v > 0 )
1515             {
1516               const int left_col = (i == 0) ? 0 : 1;
1517               tab_hline (tbl, TAL_1, left_col, n_cols - 1,
1518                          heading_rows + row_var_start + i * DESCRIPTIVE_ROWS);
1519             }
1520
1521           if ( fctr->indep_var[0])
1522             {
1523               struct string vstr;
1524               ds_init_empty (&vstr);
1525               var_append_value_name (fctr->indep_var[0],
1526                                      result->value[0], &vstr);
1527
1528               tab_text (tbl, 1,
1529                         heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1530                         TAB_LEFT,
1531                         ds_cstr (&vstr)
1532                         );
1533
1534               ds_destroy (&vstr);
1535             }
1536
1537
1538           tab_text (tbl, n_cols - 4,
1539                     heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1540                     TAB_LEFT,
1541                     _("Mean"));
1542
1543           tab_text (tbl, n_cols - 4,
1544                     heading_rows + row_var_start + 1 + i * DESCRIPTIVE_ROWS,
1545                     TAB_LEFT | TAT_PRINTF,
1546                     _("%g%% Confidence Interval for Mean"),
1547                     cmd.n_cinterval[0]);
1548
1549           tab_text (tbl, n_cols - 3,
1550                     heading_rows + row_var_start + 1 + i * DESCRIPTIVE_ROWS,
1551                     TAB_LEFT,
1552                     _("Lower Bound"));
1553
1554           tab_text (tbl, n_cols - 3,
1555                     heading_rows + row_var_start + 2 + i * DESCRIPTIVE_ROWS,
1556                     TAB_LEFT,
1557                     _("Upper Bound"));
1558
1559           tab_text (tbl, n_cols - 4,
1560                     heading_rows + row_var_start + 3 + i * DESCRIPTIVE_ROWS,
1561                     TAB_LEFT | TAT_PRINTF,
1562                     _("5%% Trimmed Mean"));
1563
1564           tab_text (tbl, n_cols - 4,
1565                     heading_rows + row_var_start + 4 + i * DESCRIPTIVE_ROWS,
1566                     TAB_LEFT,
1567                     _("Median"));
1568
1569           tab_text (tbl, n_cols - 4,
1570                     heading_rows + row_var_start + 5 + i * DESCRIPTIVE_ROWS,
1571                     TAB_LEFT,
1572                     _("Variance"));
1573
1574           tab_text (tbl, n_cols - 4,
1575                     heading_rows + row_var_start + 6 + i * DESCRIPTIVE_ROWS,
1576                     TAB_LEFT,
1577                     _("Std. Deviation"));
1578
1579           tab_text (tbl, n_cols - 4,
1580                     heading_rows + row_var_start + 7 + i * DESCRIPTIVE_ROWS,
1581                     TAB_LEFT,
1582                     _("Minimum"));
1583
1584           tab_text (tbl, n_cols - 4,
1585                     heading_rows + row_var_start + 8 + i * DESCRIPTIVE_ROWS,
1586                     TAB_LEFT,
1587                     _("Maximum"));
1588
1589           tab_text (tbl, n_cols - 4,
1590                     heading_rows + row_var_start + 9 + i * DESCRIPTIVE_ROWS,
1591                     TAB_LEFT,
1592                     _("Range"));
1593
1594           tab_text (tbl, n_cols - 4,
1595                     heading_rows + row_var_start + 10 + i * DESCRIPTIVE_ROWS,
1596                     TAB_LEFT,
1597                     _("Interquartile Range"));
1598
1599
1600           tab_text (tbl, n_cols - 4,
1601                     heading_rows + row_var_start + 11 + i * DESCRIPTIVE_ROWS,
1602                     TAB_LEFT,
1603                     _("Skewness"));
1604
1605           tab_text (tbl, n_cols - 4,
1606                     heading_rows + row_var_start + 12 + i * DESCRIPTIVE_ROWS,
1607                     TAB_LEFT,
1608                     _("Kurtosis"));
1609
1610
1611           /* Now the statistics ... */
1612
1613           tab_float (tbl, n_cols - 2,
1614                     heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1615                      TAB_CENTER,
1616                      result->metrics[v].mean,
1617                      8, 2);
1618
1619           tab_float (tbl, n_cols - 1,
1620                     heading_rows + row_var_start + i * DESCRIPTIVE_ROWS,
1621                      TAB_CENTER,
1622                      result->metrics[v].se_mean,
1623                      8, 3);
1624
1625
1626           tab_float (tbl, n_cols - 2,
1627                      heading_rows + row_var_start + 1 + i * DESCRIPTIVE_ROWS,
1628                      TAB_CENTER,
1629                      result->metrics[v].mean - t *
1630                       result->metrics[v].se_mean,
1631                      8, 3);
1632
1633           tab_float (tbl, n_cols - 2,
1634                      heading_rows + row_var_start + 2 + i * DESCRIPTIVE_ROWS,
1635                      TAB_CENTER,
1636                      result->metrics[v].mean + t *
1637                       result->metrics[v].se_mean,
1638                      8, 3);
1639
1640
1641           tab_float (tbl, n_cols - 2,
1642                      heading_rows + row_var_start + 3 + i * DESCRIPTIVE_ROWS,
1643                      TAB_CENTER,
1644                      trimmed_mean_calculate ((struct trimmed_mean *) result->metrics[v].trimmed_mean),
1645                      8, 2);
1646
1647
1648           tab_float (tbl, n_cols - 2,
1649                      heading_rows + row_var_start + 4 + i * DESCRIPTIVE_ROWS,
1650                      TAB_CENTER,
1651                      percentile_calculate (result->metrics[v].quartiles[1], percentile_algorithm),
1652                      8, 2);
1653
1654
1655           tab_float (tbl, n_cols - 2,
1656                      heading_rows + row_var_start + 5 + i * DESCRIPTIVE_ROWS,
1657                      TAB_CENTER,
1658                      result->metrics[v].variance,
1659                      8, 3);
1660
1661           tab_float (tbl, n_cols - 2,
1662                      heading_rows + row_var_start + 6 + i * DESCRIPTIVE_ROWS,
1663                      TAB_CENTER,
1664                      sqrt (result->metrics[v].variance),
1665                      8, 3);
1666
1667           tab_float (tbl, n_cols - 2,
1668                      heading_rows + row_var_start + 10 + i * DESCRIPTIVE_ROWS,
1669                      TAB_CENTER,
1670                      percentile_calculate (result->metrics[v].quartiles[2],
1671                                            percentile_algorithm) -
1672                      percentile_calculate (result->metrics[v].quartiles[0],
1673                                            percentile_algorithm),
1674                      8, 2);
1675
1676
1677           tab_float (tbl, n_cols - 2,
1678                      heading_rows + row_var_start + 11 + i * DESCRIPTIVE_ROWS,
1679                      TAB_CENTER,
1680                      result->metrics[v].skewness,
1681                      8, 3);
1682
1683           tab_float (tbl, n_cols - 2,
1684                      heading_rows + row_var_start + 12 + i * DESCRIPTIVE_ROWS,
1685                      TAB_CENTER,
1686                      result->metrics[v].kurtosis,
1687                      8, 3);
1688
1689           tab_float (tbl, n_cols - 1,
1690                      heading_rows + row_var_start + 11 + i * DESCRIPTIVE_ROWS,
1691                      TAB_CENTER,
1692                      calc_seskew (result->metrics[v].n),
1693                      8, 3);
1694
1695           tab_float (tbl, n_cols - 1,
1696                      heading_rows + row_var_start + 12 + i * DESCRIPTIVE_ROWS,
1697                      TAB_CENTER,
1698                      calc_sekurt (result->metrics[v].n),
1699                      8, 3);
1700
1701           {
1702             struct extremum *minimum, *maximum ;
1703
1704             struct ll *max_ll = ll_head (extrema_list (result->metrics[v].maxima));
1705             struct ll *min_ll = ll_head (extrema_list (result->metrics[v].minima));
1706
1707             maximum = ll_data (max_ll, struct extremum, ll);
1708             minimum = ll_data (min_ll, struct extremum, ll);
1709
1710             tab_float (tbl, n_cols - 2,
1711                        heading_rows + row_var_start + 7 + i * DESCRIPTIVE_ROWS,
1712                        TAB_CENTER,
1713                        minimum->value,
1714                        8, 3);
1715
1716             tab_float (tbl, n_cols - 2,
1717                        heading_rows + row_var_start + 8 + i * DESCRIPTIVE_ROWS,
1718                        TAB_CENTER,
1719                        maximum->value,
1720                        8, 3);
1721
1722             tab_float (tbl, n_cols - 2,
1723                        heading_rows + row_var_start + 9 + i * DESCRIPTIVE_ROWS,
1724                        TAB_CENTER,
1725                        maximum->value - minimum->value,
1726                        8, 3);
1727           }
1728         }
1729     }
1730
1731   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
1732
1733   tab_title (tbl, _("Descriptives"));
1734
1735   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE,
1736             _("Statistic"));
1737
1738   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE,
1739             _("Std. Error"));
1740
1741   tab_submit (tbl);
1742 }
1743
1744
1745
1746 static void
1747 show_extremes (const struct variable **dependent_var,
1748                int n_dep_var,
1749                const struct xfactor *fctr)
1750 {
1751   int v;
1752   int heading_columns = 3;
1753   int n_cols;
1754   const int heading_rows = 1;
1755   struct tab_table *tbl;
1756
1757   int n_rows ;
1758   n_rows = n_dep_var;
1759
1760   assert (fctr);
1761
1762   if ( fctr->indep_var[0] )
1763     {
1764       heading_columns = 4;
1765
1766       if ( fctr->indep_var[1] )
1767         {
1768           heading_columns = 5;
1769         }
1770     }
1771
1772   n_rows *= ll_count (&fctr->result_list) * cmd.st_n * 2;
1773   n_rows += heading_rows;
1774
1775   n_cols = heading_columns + 2;
1776
1777   tbl = tab_create (n_cols, n_rows, 0);
1778   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1779
1780   tab_dim (tbl, tab_natural_dimensions);
1781
1782   /* Outline the box */
1783   tab_box (tbl,
1784            TAL_2, TAL_2,
1785            -1, -1,
1786            0, 0,
1787            n_cols - 1, n_rows - 1);
1788
1789
1790   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1791   tab_hline (tbl, TAL_2, 1, n_cols - 1, heading_rows );
1792   tab_vline (tbl, TAL_1, n_cols - 1, 0, n_rows - 1);
1793
1794   if ( fctr->indep_var[0])
1795     tab_text (tbl, 1, 0, TAT_TITLE, var_to_string (fctr->indep_var[0]));
1796
1797   if ( fctr->indep_var[1])
1798     tab_text (tbl, 2, 0, TAT_TITLE, var_to_string (fctr->indep_var[1]));
1799
1800   for (v = 0 ; v < n_dep_var ; ++v )
1801     {
1802       struct ll *ll;
1803       int i = 0;
1804       const int row_var_start = v * cmd.st_n * 2 * ll_count(&fctr->result_list);
1805
1806       tab_text (tbl,
1807                 0,
1808                 heading_rows + row_var_start,
1809                 TAB_LEFT | TAT_TITLE,
1810                 var_to_string (dependent_var[v])
1811                 );
1812
1813       for (ll = ll_head (&fctr->result_list);
1814            ll != ll_null (&fctr->result_list); i++, ll = ll_next (ll))
1815         {
1816           int e ;
1817           struct ll *min_ll;
1818           struct ll *max_ll;
1819           const int row_result_start = i * cmd.st_n * 2;
1820
1821           const struct factor_result *result =
1822             ll_data (ll, struct factor_result, ll);
1823
1824           if (i > 0 || v > 0)
1825             tab_hline (tbl, TAL_1, 1, n_cols - 1,
1826                        heading_rows + row_var_start + row_result_start);
1827
1828           tab_hline (tbl, TAL_1, heading_columns - 2, n_cols - 1,
1829                      heading_rows + row_var_start + row_result_start + cmd.st_n);
1830
1831           for ( e = 1; e <= cmd.st_n; ++e )
1832             {
1833               tab_text (tbl, n_cols - 3,
1834                         heading_rows + row_var_start + row_result_start + e - 1,
1835                         TAB_RIGHT | TAT_PRINTF,
1836                         _("%d"), e);
1837
1838               tab_text (tbl, n_cols - 3,
1839                         heading_rows + row_var_start + row_result_start + cmd.st_n + e - 1,
1840                         TAB_RIGHT | TAT_PRINTF,
1841                         _("%d"), e);
1842             }
1843
1844
1845           min_ll = ll_head (extrema_list (result->metrics[v].minima));
1846           for (e = 0; e < cmd.st_n;)
1847             {
1848               struct extremum *minimum = ll_data (min_ll, struct extremum, ll);
1849               double weight = minimum->weight;
1850
1851               while (weight-- > 0 && e < cmd.st_n)
1852                 {
1853                   tab_float (tbl, n_cols - 1,
1854                              heading_rows + row_var_start + row_result_start + cmd.st_n + e,
1855                              TAB_RIGHT,
1856                              minimum->value,
1857                              8, 2);
1858
1859
1860                   tab_float (tbl, n_cols - 2,
1861                              heading_rows + row_var_start + row_result_start + cmd.st_n + e,
1862                              TAB_RIGHT,
1863                              minimum->location,
1864                              8, 0);
1865                   ++e;
1866                 }
1867
1868               min_ll = ll_next (min_ll);
1869             }
1870
1871
1872           max_ll = ll_head (extrema_list (result->metrics[v].maxima));
1873           for (e = 0; e < cmd.st_n;)
1874             {
1875               struct extremum *maximum = ll_data (max_ll, struct extremum, ll);
1876               double weight = maximum->weight;
1877
1878               while (weight-- > 0 && e < cmd.st_n)
1879                 {
1880                   tab_float (tbl, n_cols - 1,
1881                              heading_rows + row_var_start + row_result_start + e,
1882                              TAB_RIGHT,
1883                              maximum->value,
1884                              8, 2);
1885
1886
1887                   tab_float (tbl, n_cols - 2,
1888                              heading_rows + row_var_start + row_result_start + e,
1889                              TAB_RIGHT,
1890                              maximum->location,
1891                              8, 0);
1892                   ++e;
1893                 }
1894
1895               max_ll = ll_next (max_ll);
1896             }
1897
1898
1899           if ( fctr->indep_var[0])
1900             {
1901               struct string vstr;
1902               ds_init_empty (&vstr);
1903               var_append_value_name (fctr->indep_var[0],
1904                                      result->value[0], &vstr);
1905
1906               tab_text (tbl, 1,
1907                         heading_rows + row_var_start + row_result_start,
1908                         TAB_LEFT,
1909                         ds_cstr (&vstr)
1910                         );
1911
1912               ds_destroy (&vstr);
1913             }
1914
1915
1916           tab_text (tbl, n_cols - 4,
1917                     heading_rows + row_var_start + row_result_start,
1918                     TAB_RIGHT,
1919                     _("Highest"));
1920
1921           tab_text (tbl, n_cols - 4,
1922                     heading_rows + row_var_start + row_result_start + cmd.st_n,
1923                     TAB_RIGHT,
1924                     _("Lowest"));
1925         }
1926     }
1927
1928   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
1929
1930
1931   tab_title (tbl, _("Extreme Values"));
1932
1933
1934   tab_text (tbl, n_cols - 2, 0, TAB_CENTER | TAT_TITLE,
1935             _("Case Number"));
1936
1937
1938   tab_text (tbl, n_cols - 1, 0, TAB_CENTER | TAT_TITLE,
1939             _("Value"));
1940
1941   tab_submit (tbl);
1942 }
1943
1944 #define PERCENTILE_ROWS 2
1945
1946 static void
1947 show_percentiles (const struct variable **dependent_var,
1948                   int n_dep_var,
1949                   const struct xfactor *fctr)
1950 {
1951   int i;
1952   int v;
1953   int heading_columns = 2;
1954   int n_cols;
1955   const int n_percentiles = subc_list_double_count (&percentile_list);
1956   const int heading_rows = 2;
1957   struct tab_table *tbl;
1958
1959   int n_rows ;
1960   n_rows = n_dep_var;
1961
1962   assert (fctr);
1963
1964   if ( fctr->indep_var[0] )
1965     {
1966       heading_columns = 3;
1967
1968       if ( fctr->indep_var[1] )
1969         {
1970           heading_columns = 4;
1971         }
1972     }
1973
1974   n_rows *= ll_count (&fctr->result_list) * PERCENTILE_ROWS;
1975   n_rows += heading_rows;
1976
1977   n_cols = heading_columns + n_percentiles;
1978
1979   tbl = tab_create (n_cols, n_rows, 0);
1980   tab_headers (tbl, heading_columns, 0, heading_rows, 0);
1981
1982   tab_dim (tbl, tab_natural_dimensions);
1983
1984   /* Outline the box */
1985   tab_box (tbl,
1986            TAL_2, TAL_2,
1987            -1, -1,
1988            0, 0,
1989            n_cols - 1, n_rows - 1);
1990
1991
1992   tab_hline (tbl, TAL_2, 0, n_cols - 1, heading_rows );
1993   tab_hline (tbl, TAL_2, 1, n_cols - 1, heading_rows );
1994
1995   if ( fctr->indep_var[0])
1996     tab_text (tbl, 1, 1, TAT_TITLE, var_to_string (fctr->indep_var[0]));
1997
1998   if ( fctr->indep_var[1])
1999     tab_text (tbl, 2, 1, TAT_TITLE, var_to_string (fctr->indep_var[1]));
2000
2001   for (v = 0 ; v < n_dep_var ; ++v )
2002     {
2003       double hinges[3];
2004       struct ll *ll;
2005       int i = 0;
2006
2007       const int row_var_start =
2008         v * PERCENTILE_ROWS * ll_count(&fctr->result_list);
2009
2010       tab_text (tbl,
2011                 0,
2012                 heading_rows + row_var_start,
2013                 TAB_LEFT | TAT_TITLE,
2014                 var_to_string (dependent_var[v])
2015                 );
2016
2017       for (ll = ll_head (&fctr->result_list);
2018            ll != ll_null (&fctr->result_list); i++, ll = ll_next (ll))
2019         {
2020           int j;
2021           const struct factor_result *result =
2022             ll_data (ll, struct factor_result, ll);
2023
2024           if ( i > 0 || v > 0 )
2025             {
2026               const int left_col = (i == 0) ? 0 : 1;
2027               tab_hline (tbl, TAL_1, left_col, n_cols - 1,
2028                          heading_rows + row_var_start + i * PERCENTILE_ROWS);
2029             }
2030
2031           if ( fctr->indep_var[0])
2032             {
2033               struct string vstr;
2034               ds_init_empty (&vstr);
2035               var_append_value_name (fctr->indep_var[0],
2036                                      result->value[0], &vstr);
2037
2038               tab_text (tbl, 1,
2039                         heading_rows + row_var_start + i * PERCENTILE_ROWS,
2040                         TAB_LEFT,
2041                         ds_cstr (&vstr)
2042                         );
2043
2044               ds_destroy (&vstr);
2045             }
2046
2047
2048           tab_text (tbl, n_cols - n_percentiles - 1,
2049                     heading_rows + row_var_start + i * PERCENTILE_ROWS,
2050                     TAB_LEFT,
2051                     ptile_alg_desc [percentile_algorithm]);
2052
2053
2054           tab_text (tbl, n_cols - n_percentiles - 1,
2055                     heading_rows + row_var_start + 1 + i * PERCENTILE_ROWS,
2056                     TAB_LEFT,
2057                     _("Tukey's Hinges"));
2058
2059
2060           tab_vline (tbl, TAL_1, n_cols - n_percentiles -1, heading_rows, n_rows - 1);
2061
2062           tukey_hinges_calculate ((struct tukey_hinges *) result->metrics[v].tukey_hinges,
2063                                   hinges);
2064
2065           for (j = 0; j < n_percentiles; ++j)
2066             {
2067               double hinge = SYSMIS;
2068               tab_float (tbl, n_cols - n_percentiles + j,
2069                          heading_rows + row_var_start + i * PERCENTILE_ROWS,
2070                          TAB_CENTER,
2071                          percentile_calculate (result->metrics[v].ptl[j],
2072                                                percentile_algorithm),
2073                          8, 2
2074                          );
2075
2076               if ( result->metrics[v].ptl[j]->ptile == 0.5)
2077                 hinge = hinges[1];
2078               else if ( result->metrics[v].ptl[j]->ptile == 0.25)
2079                 hinge = hinges[0];
2080               else if ( result->metrics[v].ptl[j]->ptile == 0.75)
2081                 hinge = hinges[2];
2082
2083               if ( hinge != SYSMIS)
2084                 tab_float (tbl, n_cols - n_percentiles + j,
2085                            heading_rows + row_var_start + 1 + i * PERCENTILE_ROWS,
2086                            TAB_CENTER,
2087                            hinge,
2088                            8, 2
2089                            );
2090
2091             }
2092         }
2093     }
2094
2095   tab_vline (tbl, TAL_2, heading_columns, 0, n_rows - 1);
2096
2097   tab_title (tbl, _("Percentiles"));
2098
2099
2100   for (i = 0 ; i < n_percentiles; ++i )
2101     {
2102       tab_text (tbl, n_cols - n_percentiles + i, 1,
2103                 TAB_CENTER | TAT_TITLE | TAT_PRINTF,
2104                 _("%g"),
2105                 subc_list_double_at (&percentile_list, i)
2106                 );
2107
2108
2109     }
2110
2111   tab_joint_text (tbl,
2112                   n_cols - n_percentiles, 0,
2113                   n_cols - 1, 0,
2114                   TAB_CENTER | TAT_TITLE,
2115                   _("Percentiles"));
2116
2117   /* Vertical lines for the data only */
2118   tab_box (tbl,
2119            -1, -1,
2120            -1, TAL_1,
2121            n_cols - n_percentiles, 1,
2122            n_cols - 1, n_rows - 1);
2123
2124   tab_hline (tbl, TAL_1, n_cols - n_percentiles, n_cols - 1, 1);
2125
2126
2127   tab_submit (tbl);
2128 }
2129
2130
2131 static void
2132 factor_to_string_concise (const struct xfactor *fctr,
2133                           const struct factor_result *result,
2134                           struct string *str
2135                           )
2136 {
2137   if (fctr->indep_var[0])
2138     {
2139       var_append_value_name (fctr->indep_var[0], result->value[0], str);
2140
2141       if ( fctr->indep_var[1] )
2142         {
2143           ds_put_cstr (str, ",");
2144
2145           var_append_value_name (fctr->indep_var[1], result->value[1], str);
2146
2147           ds_put_cstr (str, ")");
2148         }
2149     }
2150 }
2151
2152
2153 static void
2154 factor_to_string (const struct xfactor *fctr,
2155                   const struct factor_result *result,
2156                   struct string *str
2157                   )
2158 {
2159   if (fctr->indep_var[0])
2160     {
2161       ds_put_format (str, "(%s = ", var_get_name (fctr->indep_var[0]));
2162
2163       var_append_value_name (fctr->indep_var[0], result->value[0], str);
2164
2165       if ( fctr->indep_var[1] )
2166         {
2167           ds_put_cstr (str, ",");
2168           ds_put_format (str, "%s = ", var_get_name (fctr->indep_var[1]));
2169
2170           var_append_value_name (fctr->indep_var[1], result->value[1], str);
2171         }
2172       ds_put_cstr (str, ")");
2173     }
2174 }
2175
2176
2177
2178
2179 /*
2180   Local Variables:
2181   mode: c
2182   End:
2183 */