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