CROSSTABS: Handle case where all cases in a crosstabulation are missing.
[pspp] / src / language / stats / descriptives.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 1997-9, 2000, 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 <limits.h>
20 #include <math.h>
21 #include <stdlib.h>
22
23 #include <data/casegrouper.h>
24 #include <data/casereader.h>
25 #include <data/dictionary.h>
26 #include <data/procedure.h>
27 #include <data/transformations.h>
28 #include <data/variable.h>
29 #include <language/command.h>
30 #include <language/dictionary/split-file.h>
31 #include <language/lexer/lexer.h>
32 #include <language/lexer/variable-parser.h>
33 #include <libpspp/array.h>
34 #include <libpspp/compiler.h>
35 #include <libpspp/message.h>
36 #include <libpspp/assertion.h>
37 #include <math/moments.h>
38 #include <output/tab.h>
39
40 #include "xalloc.h"
41
42 #include "gettext.h"
43 #define _(msgid) gettext (msgid)
44 #define N_(msgid) msgid
45
46 /* DESCRIPTIVES private data. */
47
48 struct dsc_proc;
49
50 /* Handling of missing values. */
51 enum dsc_missing_type
52   {
53     DSC_VARIABLE,       /* Handle missing values on a per-variable basis. */
54     DSC_LISTWISE        /* Discard entire case if any variable is missing. */
55   };
56
57 /* Describes properties of a distribution for the purpose of
58    calculating a Z-score. */
59 struct dsc_z_score
60   {
61     const struct variable *src_var;   /* Variable on which z-score is based. */
62     struct variable *z_var;     /* New z-score variable. */
63     double mean;                /* Distribution mean. */
64     double std_dev;             /* Distribution standard deviation. */
65   };
66
67 /* DESCRIPTIVES transformation (for calculating Z-scores). */
68 struct dsc_trns
69   {
70     struct dsc_z_score *z_scores; /* Array of Z-scores. */
71     int z_score_cnt;            /* Number of Z-scores. */
72     const struct variable **vars;     /* Variables for listwise missing checks. */
73     size_t var_cnt;             /* Number of variables. */
74     enum dsc_missing_type missing_type; /* Treatment of missing values. */
75     enum mv_class exclude;      /* Classes of missing values to exclude. */
76   };
77
78 /* Statistics.  Used as bit indexes, so must be 32 or fewer. */
79 enum dsc_statistic
80   {
81     DSC_MEAN = 0, DSC_SEMEAN, DSC_STDDEV, DSC_VARIANCE, DSC_KURTOSIS,
82     DSC_SEKURT, DSC_SKEWNESS, DSC_SESKEW, DSC_RANGE, DSC_MIN,
83     DSC_MAX, DSC_SUM, DSC_N_STATS,
84
85     /* Only valid as sort criteria. */
86     DSC_NAME = -2,              /* Sort by name. */
87     DSC_NONE = -1               /* Unsorted. */
88   };
89
90 /* Describes one statistic. */
91 struct dsc_statistic_info
92   {
93     const char *identifier;     /* Identifier. */
94     const char *name;           /* Full name. */
95     enum moment moment;         /* Highest moment needed to calculate. */
96   };
97
98 /* Table of statistics, indexed by DSC_*. */
99 static const struct dsc_statistic_info dsc_info[DSC_N_STATS] =
100   {
101     {"MEAN", N_("Mean"), MOMENT_MEAN},
102     {"SEMEAN", N_("S.E. Mean"), MOMENT_VARIANCE},
103     {"STDDEV", N_("Std Dev"), MOMENT_VARIANCE},
104     {"VARIANCE", N_("Variance"), MOMENT_VARIANCE},
105     {"KURTOSIS", N_("Kurtosis"), MOMENT_KURTOSIS},
106     {"SEKURTOSIS", N_("S.E. Kurt"), MOMENT_NONE},
107     {"SKEWNESS", N_("Skewness"), MOMENT_SKEWNESS},
108     {"SESKEWNESS", N_("S.E. Skew"), MOMENT_NONE},
109     {"RANGE", N_("Range"), MOMENT_NONE},
110     {"MINIMUM", N_("Minimum"), MOMENT_NONE},
111     {"MAXIMUM", N_("Maximum"), MOMENT_NONE},
112     {"SUM", N_("Sum"), MOMENT_MEAN},
113   };
114
115 /* Statistics calculated by default if none are explicitly
116    requested. */
117 #define DEFAULT_STATS                                                   \
118         ((1ul << DSC_MEAN) | (1ul << DSC_STDDEV) | (1ul << DSC_MIN)     \
119          | (1ul << DSC_MAX))
120
121 /* A variable specified on DESCRIPTIVES. */
122 struct dsc_var
123   {
124     const struct variable *v;         /* Variable to calculate on. */
125     char z_name[VAR_NAME_LEN + 1]; /* Name for z-score variable. */
126     double valid, missing;      /* Valid, missing counts. */
127     struct moments *moments;    /* Moments. */
128     double min, max;            /* Maximum and mimimum values. */
129     double stats[DSC_N_STATS];  /* All the stats' values. */
130   };
131
132 /* Output format. */
133 enum dsc_format
134   {
135     DSC_LINE,           /* Abbreviated format. */
136     DSC_SERIAL          /* Long format. */
137   };
138
139 /* A DESCRIPTIVES procedure. */
140 struct dsc_proc
141   {
142     /* Per-variable info. */
143     struct dsc_var *vars;       /* Variables. */
144     size_t var_cnt;             /* Number of variables. */
145
146     /* User options. */
147     enum dsc_missing_type missing_type; /* Treatment of missing values. */
148     enum mv_class exclude;      /* Classes of missing values to exclude. */
149     int show_var_labels;        /* Nonzero to show variable labels. */
150     int show_index;             /* Nonzero to show variable index. */
151     enum dsc_format format;     /* Output format. */
152
153     /* Accumulated results. */
154     double missing_listwise;    /* Sum of weights of cases missing listwise. */
155     double valid;               /* Sum of weights of valid cases. */
156     bool bad_warn;               /* Warn if bad weight found. */
157     enum dsc_statistic sort_by_stat; /* Statistic to sort by; -1: name. */
158     int sort_ascending;         /* !0: ascending order; 0: descending. */
159     unsigned long show_stats;   /* Statistics to display. */
160     unsigned long calc_stats;   /* Statistics to calculate. */
161     enum moment max_moment;     /* Highest moment needed for stats. */
162   };
163
164 /* Parsing. */
165 static enum dsc_statistic match_statistic (struct lexer *);
166 static void free_dsc_proc (struct dsc_proc *);
167
168 /* Z-score functions. */
169 static bool try_name (const struct dictionary *dict,
170                       struct dsc_proc *dsc, const char *name);
171 static bool generate_z_varname (const struct dictionary *dict,
172                                 struct dsc_proc *dsc, char *z_name,
173                                 const char *name, int *z_cnt);
174 static void dump_z_table (struct dsc_proc *);
175 static void setup_z_trns (struct dsc_proc *, struct dataset *);
176
177 /* Procedure execution functions. */
178 static void calc_descriptives (struct dsc_proc *, struct casereader *,
179                                struct dataset *);
180 static void display (struct dsc_proc *dsc);
181 \f
182 /* Parser and outline. */
183
184 /* Handles DESCRIPTIVES. */
185 int
186 cmd_descriptives (struct lexer *lexer, struct dataset *ds)
187 {
188   struct dictionary *dict = dataset_dict (ds);
189   struct dsc_proc *dsc;
190   const struct variable **vars = NULL;
191   size_t var_cnt = 0;
192   int save_z_scores = 0;
193   int z_cnt = 0;
194   size_t i;
195   bool ok;
196
197   struct casegrouper *grouper;
198   struct casereader *group;
199
200   /* Create and initialize dsc. */
201   dsc = xmalloc (sizeof *dsc);
202   dsc->vars = NULL;
203   dsc->var_cnt = 0;
204   dsc->missing_type = DSC_VARIABLE;
205   dsc->exclude = MV_ANY;
206   dsc->show_var_labels = 1;
207   dsc->show_index = 0;
208   dsc->format = DSC_LINE;
209   dsc->missing_listwise = 0.;
210   dsc->valid = 0.;
211   dsc->bad_warn = 1;
212   dsc->sort_by_stat = DSC_NONE;
213   dsc->sort_ascending = 1;
214   dsc->show_stats = dsc->calc_stats = DEFAULT_STATS;
215
216   /* Parse DESCRIPTIVES. */
217   while (lex_token (lexer) != '.')
218     {
219       if (lex_match_id (lexer, "MISSING"))
220         {
221           lex_match (lexer, '=');
222           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
223             {
224               if (lex_match_id (lexer, "VARIABLE"))
225                 dsc->missing_type = DSC_VARIABLE;
226               else if (lex_match_id (lexer, "LISTWISE"))
227                 dsc->missing_type = DSC_LISTWISE;
228               else if (lex_match_id (lexer, "INCLUDE"))
229                 dsc->exclude = MV_SYSTEM;
230               else
231                 {
232                   lex_error (lexer, NULL);
233                   goto error;
234                 }
235               lex_match (lexer, ',');
236             }
237         }
238       else if (lex_match_id (lexer, "SAVE"))
239         save_z_scores = 1;
240       else if (lex_match_id (lexer, "FORMAT"))
241         {
242           lex_match (lexer, '=');
243           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
244             {
245               if (lex_match_id (lexer, "LABELS"))
246                 dsc->show_var_labels = 1;
247               else if (lex_match_id (lexer, "NOLABELS"))
248                 dsc->show_var_labels = 0;
249               else if (lex_match_id (lexer, "INDEX"))
250                 dsc->show_index = 1;
251               else if (lex_match_id (lexer, "NOINDEX"))
252                 dsc->show_index = 0;
253               else if (lex_match_id (lexer, "LINE"))
254                 dsc->format = DSC_LINE;
255               else if (lex_match_id (lexer, "SERIAL"))
256                 dsc->format = DSC_SERIAL;
257               else
258                 {
259                   lex_error (lexer, NULL);
260                   goto error;
261                 }
262               lex_match (lexer, ',');
263             }
264         }
265       else if (lex_match_id (lexer, "STATISTICS"))
266         {
267           lex_match (lexer, '=');
268           dsc->show_stats = 0;
269           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
270             {
271               if (lex_match (lexer, T_ALL))
272                 dsc->show_stats |= (1ul << DSC_N_STATS) - 1;
273               else if (lex_match_id (lexer, "DEFAULT"))
274                 dsc->show_stats |= DEFAULT_STATS;
275               else
276                 dsc->show_stats |= 1ul << (match_statistic (lexer));
277               lex_match (lexer, ',');
278             }
279           if (dsc->show_stats == 0)
280             dsc->show_stats = DEFAULT_STATS;
281         }
282       else if (lex_match_id (lexer, "SORT"))
283         {
284           lex_match (lexer, '=');
285           if (lex_match_id (lexer, "NAME"))
286             dsc->sort_by_stat = DSC_NAME;
287           else
288             {
289               dsc->sort_by_stat = match_statistic (lexer);
290               if (dsc->sort_by_stat == DSC_NONE )
291                 dsc->sort_by_stat = DSC_MEAN;
292             }
293           if (lex_match (lexer, '('))
294             {
295               if (lex_match_id (lexer, "A"))
296                 dsc->sort_ascending = 1;
297               else if (lex_match_id (lexer, "D"))
298                 dsc->sort_ascending = 0;
299               else
300                 lex_error (lexer, NULL);
301               lex_force_match (lexer, ')');
302             }
303         }
304       else if (var_cnt == 0)
305         {
306           if (lex_look_ahead (lexer) == '=')
307             {
308               lex_match_id (lexer, "VARIABLES");
309               lex_match (lexer, '=');
310             }
311
312           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
313             {
314               int i;
315
316               if (!parse_variables_const (lexer, dict, &vars, &var_cnt,
317                                     PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
318                 goto error;
319
320               dsc->vars = xnrealloc ((void *)dsc->vars, var_cnt, sizeof *dsc->vars);
321               for (i = dsc->var_cnt; i < var_cnt; i++)
322                 {
323                   struct dsc_var *dv = &dsc->vars[i];
324                   dv->v = vars[i];
325                   dv->z_name[0] = '\0';
326                   dv->moments = NULL;
327                 }
328               dsc->var_cnt = var_cnt;
329
330               if (lex_match (lexer, '('))
331                 {
332                   if (lex_token (lexer) != T_ID)
333                     {
334                       lex_error (lexer, NULL);
335                       goto error;
336                     }
337                   if (try_name (dict, dsc, lex_tokid (lexer)))
338                     {
339                       strcpy (dsc->vars[dsc->var_cnt - 1].z_name, lex_tokid (lexer));
340                       z_cnt++;
341                     }
342                   else
343                     msg (SE, _("Z-score variable name %s would be"
344                                " a duplicate variable name."), lex_tokid (lexer));
345                   lex_get (lexer);
346                   if (!lex_force_match (lexer, ')'))
347                     goto error;
348                 }
349             }
350         }
351       else
352         {
353           lex_error (lexer, NULL);
354           goto error;
355         }
356
357       lex_match (lexer, '/');
358     }
359   if (var_cnt == 0)
360     {
361       msg (SE, _("No variables specified."));
362       goto error;
363     }
364
365   /* Construct z-score varnames, show translation table. */
366   if (z_cnt || save_z_scores)
367     {
368       if (save_z_scores)
369         {
370           int gen_cnt = 0;
371
372           for (i = 0; i < dsc->var_cnt; i++)
373             if (dsc->vars[i].z_name[0] == 0)
374               {
375                 if (!generate_z_varname (dict, dsc, dsc->vars[i].z_name,
376                                          var_get_name (dsc->vars[i].v),
377                                          &gen_cnt))
378                   goto error;
379                 z_cnt++;
380               }
381         }
382       dump_z_table (dsc);
383     }
384
385   /* Figure out statistics to display. */
386   if (dsc->show_stats & (1ul << DSC_SKEWNESS))
387     dsc->show_stats |= 1ul << DSC_SESKEW;
388   if (dsc->show_stats & (1ul << DSC_KURTOSIS))
389     dsc->show_stats |= 1ul << DSC_SEKURT;
390
391   /* Figure out which statistics to calculate. */
392   dsc->calc_stats = dsc->show_stats;
393   if (z_cnt > 0)
394     dsc->calc_stats |= (1ul << DSC_MEAN) | (1ul << DSC_STDDEV);
395   if (dsc->sort_by_stat >= 0)
396     dsc->calc_stats |= 1ul << dsc->sort_by_stat;
397   if (dsc->show_stats & (1ul << DSC_SESKEW))
398     dsc->calc_stats |= 1ul << DSC_SKEWNESS;
399   if (dsc->show_stats & (1ul << DSC_SEKURT))
400     dsc->calc_stats |= 1ul << DSC_KURTOSIS;
401
402   /* Figure out maximum moment needed and allocate moments for
403      the variables. */
404   dsc->max_moment = MOMENT_NONE;
405   for (i = 0; i < DSC_N_STATS; i++)
406     if (dsc->calc_stats & (1ul << i) && dsc_info[i].moment > dsc->max_moment)
407       dsc->max_moment = dsc_info[i].moment;
408   if (dsc->max_moment != MOMENT_NONE)
409     for (i = 0; i < dsc->var_cnt; i++)
410       dsc->vars[i].moments = moments_create (dsc->max_moment);
411
412   /* Data pass. */
413   grouper = casegrouper_create_splits (proc_open (ds), dict);
414   while (casegrouper_get_next_group (grouper, &group))
415     calc_descriptives (dsc, group, ds);
416   ok = casegrouper_destroy (grouper);
417   ok = proc_commit (ds) && ok;
418
419   /* Z-scoring! */
420   if (ok && z_cnt)
421     setup_z_trns (dsc, ds);
422
423   /* Done. */
424   free (vars);
425   free_dsc_proc (dsc);
426   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
427
428  error:
429   free (vars);
430   free_dsc_proc (dsc);
431   return CMD_FAILURE;
432 }
433
434 /* Returns the statistic named by the current token and skips past the token.
435    Returns DSC_NONE if no statistic is given (e.g., subcommand with no
436    specifiers). Emits an error if the current token ID does not name a
437    statistic. */
438 static enum dsc_statistic
439 match_statistic (struct lexer *lexer)
440 {
441   if (lex_token (lexer) == T_ID)
442     {
443       enum dsc_statistic stat;
444
445       for (stat = 0; stat < DSC_N_STATS; stat++)
446         if (lex_match_id (lexer, dsc_info[stat].identifier))
447           return stat;
448
449       lex_get (lexer);
450       lex_error (lexer, _("expecting statistic name: reverting to default"));
451     }
452
453   return DSC_NONE;
454 }
455
456 /* Frees DSC. */
457 static void
458 free_dsc_proc (struct dsc_proc *dsc)
459 {
460   size_t i;
461
462   if (dsc == NULL)
463     return;
464
465   for (i = 0; i < dsc->var_cnt; i++)
466     moments_destroy (dsc->vars[i].moments);
467   free (dsc->vars);
468   free (dsc);
469 }
470 \f
471 /* Z scores. */
472
473 /* Returns false if NAME is a duplicate of any existing variable name or
474    of any previously-declared z-var name; otherwise returns true. */
475 static bool
476 try_name (const struct dictionary *dict, struct dsc_proc *dsc,
477           const char *name)
478 {
479   size_t i;
480
481   if (dict_lookup_var (dict, name) != NULL)
482     return false;
483   for (i = 0; i < dsc->var_cnt; i++)
484     if (!strcasecmp (dsc->vars[i].z_name, name))
485       return false;
486   return true;
487 }
488
489 /* Generates a name for a Z-score variable based on a variable
490    named VAR_NAME, given that *Z_CNT generated variable names are
491    known to already exist.  If successful, returns true and
492    copies the new name into Z_NAME.  On failure, returns false. */
493 static bool
494 generate_z_varname (const struct dictionary *dict, struct dsc_proc *dsc, char *z_name,
495                     const char *var_name, int *z_cnt)
496 {
497   char name[VAR_NAME_LEN + 1];
498
499   /* Try a name based on the original variable name. */
500   name[0] = 'Z';
501   str_copy_trunc (name + 1, sizeof name - 1, var_name);
502   if (try_name (dict, dsc, name))
503     {
504       strcpy (z_name, name);
505       return true;
506     }
507
508   /* Generate a synthetic name. */
509   for (;;)
510     {
511       (*z_cnt)++;
512
513       if (*z_cnt <= 99)
514         sprintf (name, "ZSC%03d", *z_cnt);
515       else if (*z_cnt <= 108)
516         sprintf (name, "STDZ%02d", *z_cnt - 99);
517       else if (*z_cnt <= 117)
518         sprintf (name, "ZZZZ%02d", *z_cnt - 108);
519       else if (*z_cnt <= 126)
520         sprintf (name, "ZQZQ%02d", *z_cnt - 117);
521       else
522         {
523           msg (SE, _("Ran out of generic names for Z-score variables.  "
524                      "There are only 126 generic names: ZSC001-ZSC0999, "
525                      "STDZ01-STDZ09, ZZZZ01-ZZZZ09, ZQZQ01-ZQZQ09."));
526           return false;
527         }
528
529       if (try_name (dict, dsc, name))
530         {
531           strcpy (z_name, name);
532           return true;
533         }
534     }
535   NOT_REACHED();
536 }
537
538 /* Outputs a table describing the mapping between source
539    variables and Z-score variables. */
540 static void
541 dump_z_table (struct dsc_proc *dsc)
542 {
543   size_t cnt = 0;
544   struct tab_table *t;
545
546   {
547     size_t i;
548
549     for (i = 0; i < dsc->var_cnt; i++)
550       if (dsc->vars[i].z_name[0] != '\0')
551         cnt++;
552   }
553
554   t = tab_create (2, cnt + 1);
555   tab_title (t, _("Mapping of variables to corresponding Z-scores."));
556   tab_headers (t, 0, 0, 1, 0);
557   tab_box (t, TAL_1, TAL_1, TAL_0, TAL_1, 0, 0, 1, cnt);
558   tab_hline (t, TAL_2, 0, 1, 1);
559   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
560   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Target"));
561
562   {
563     size_t i, y;
564
565     for (i = 0, y = 1; i < dsc->var_cnt; i++)
566       if (dsc->vars[i].z_name[0] != '\0')
567         {
568           tab_text (t, 0, y, TAB_LEFT, var_get_name (dsc->vars[i].v));
569           tab_text (t, 1, y++, TAB_LEFT, dsc->vars[i].z_name);
570         }
571   }
572
573   tab_submit (t);
574 }
575
576 /* Transformation function to calculate Z-scores. Will return SYSMIS if any of
577    the following are true: 1) mean or standard deviation is SYSMIS 2) score is
578    SYSMIS 3) score is user missing and they were not included in the original
579    analyis. 4) any of the variables in the original analysis were missing
580    (either system or user-missing values that weren't included).
581 */
582 static int
583 descriptives_trns_proc (void *trns_, struct ccase **c,
584                         casenumber case_idx UNUSED)
585 {
586   struct dsc_trns *t = trns_;
587   struct dsc_z_score *z;
588   const struct variable **vars;
589   int all_sysmis = 0;
590
591   if (t->missing_type == DSC_LISTWISE)
592     {
593       assert(t->vars);
594       for (vars = t->vars; vars < t->vars + t->var_cnt; vars++)
595         {
596           double score = case_num (*c, *vars);
597           if (var_is_num_missing (*vars, score, t->exclude))
598             {
599               all_sysmis = 1;
600               break;
601             }
602         }
603     }
604
605   *c = case_unshare (*c);
606   for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
607     {
608       double input = case_num (*c, z->src_var);
609       double *output = &case_data_rw (*c, z->z_var)->f;
610
611       if (z->mean == SYSMIS || z->std_dev == SYSMIS || all_sysmis
612           || var_is_num_missing (z->src_var, input, t->exclude))
613         *output = SYSMIS;
614       else
615         *output = (input - z->mean) / z->std_dev;
616     }
617   return TRNS_CONTINUE;
618 }
619
620 /* Frees a descriptives_trns struct. */
621 static bool
622 descriptives_trns_free (void *trns_)
623 {
624   struct dsc_trns *t = trns_;
625
626   free (t->z_scores);
627   assert((t->missing_type != DSC_LISTWISE) ^ (t->vars != NULL));
628   free (t->vars);
629   return true;
630 }
631
632 /* Sets up a transformation to calculate Z scores. */
633 static void
634 setup_z_trns (struct dsc_proc *dsc, struct dataset *ds)
635 {
636   struct dsc_trns *t;
637   size_t cnt, i;
638
639   for (cnt = i = 0; i < dsc->var_cnt; i++)
640     if (dsc->vars[i].z_name[0] != '\0')
641       cnt++;
642
643   t = xmalloc (sizeof *t);
644   t->z_scores = xnmalloc (cnt, sizeof *t->z_scores);
645   t->z_score_cnt = cnt;
646   t->missing_type = dsc->missing_type;
647   t->exclude = dsc->exclude;
648   if ( t->missing_type == DSC_LISTWISE )
649     {
650       t->var_cnt = dsc->var_cnt;
651       t->vars = xnmalloc (t->var_cnt, sizeof *t->vars);
652       for (i = 0; i < t->var_cnt; i++)
653         t->vars[i] = dsc->vars[i].v;
654     }
655   else
656     {
657       t->var_cnt = 0;
658       t->vars = NULL;
659     }
660
661   for (cnt = i = 0; i < dsc->var_cnt; i++)
662     {
663       struct dsc_var *dv = &dsc->vars[i];
664       if (dv->z_name[0] != '\0')
665         {
666           struct dsc_z_score *z;
667           struct variable *dst_var;
668
669           dst_var = dict_create_var_assert (dataset_dict (ds), dv->z_name, 0);
670           var_set_label (dst_var, xasprintf (_("Z-score of %s"),
671                                              var_to_string (dv->v)));
672
673           z = &t->z_scores[cnt++];
674           z->src_var = dv->v;
675           z->z_var = dst_var;
676           z->mean = dv->stats[DSC_MEAN];
677           z->std_dev = dv->stats[DSC_STDDEV];
678         }
679     }
680
681   add_transformation (ds,
682                       descriptives_trns_proc, descriptives_trns_free, t);
683 }
684 \f
685 /* Statistical calculation. */
686
687 static bool listwise_missing (struct dsc_proc *dsc, const struct ccase *c);
688
689 /* Calculates and displays descriptive statistics for the cases
690    in CF. */
691 static void
692 calc_descriptives (struct dsc_proc *dsc, struct casereader *group,
693                    struct dataset *ds)
694 {
695   struct casereader *pass1, *pass2;
696   struct ccase *c;
697   size_t i;
698
699   c = casereader_peek (group, 0);
700   if (c == NULL)
701     {
702       casereader_destroy (group);
703       return;
704     }
705   output_split_file_values (ds, c);
706   case_unref (c);
707
708   group = casereader_create_filter_weight (group, dataset_dict (ds),
709                                            NULL, NULL);
710
711   pass1 = group;
712   pass2 = dsc->max_moment <= MOMENT_MEAN ? NULL : casereader_clone (pass1);
713
714   for (i = 0; i < dsc->var_cnt; i++)
715     {
716       struct dsc_var *dv = &dsc->vars[i];
717
718       dv->valid = dv->missing = 0.0;
719       if (dv->moments != NULL)
720         moments_clear (dv->moments);
721       dv->min = DBL_MAX;
722       dv->max = -DBL_MAX;
723     }
724   dsc->missing_listwise = 0.;
725   dsc->valid = 0.;
726
727   /* First pass to handle most of the work. */
728   for (; (c = casereader_read (pass1)) != NULL; case_unref (c))
729     {
730       double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
731
732       /* Check for missing values. */
733       if (listwise_missing (dsc, c))
734         {
735           dsc->missing_listwise += weight;
736           if (dsc->missing_type == DSC_LISTWISE)
737             continue;
738         }
739       dsc->valid += weight;
740
741       for (i = 0; i < dsc->var_cnt; i++)
742         {
743           struct dsc_var *dv = &dsc->vars[i];
744           double x = case_num (c, dv->v);
745
746           if (var_is_num_missing (dv->v, x, dsc->exclude))
747             {
748               dv->missing += weight;
749               continue;
750             }
751
752           if (dv->moments != NULL)
753             moments_pass_one (dv->moments, x, weight);
754
755           if (x < dv->min)
756             dv->min = x;
757           if (x > dv->max)
758             dv->max = x;
759         }
760     }
761   if (!casereader_destroy (pass1))
762     {
763       casereader_destroy (pass2);
764       return;
765     }
766
767   /* Second pass for higher-order moments. */
768   if (dsc->max_moment > MOMENT_MEAN)
769     {
770       for (; (c = casereader_read (pass2)) != NULL; case_unref (c))
771         {
772           double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
773
774           /* Check for missing values. */
775           if (dsc->missing_type == DSC_LISTWISE && listwise_missing (dsc, c))
776             continue;
777
778           for (i = 0; i < dsc->var_cnt; i++)
779             {
780               struct dsc_var *dv = &dsc->vars[i];
781               double x = case_num (c, dv->v);
782
783               if (var_is_num_missing (dv->v, x, dsc->exclude))
784                 continue;
785
786               if (dv->moments != NULL)
787                 moments_pass_two (dv->moments, x, weight);
788             }
789         }
790       if (!casereader_destroy (pass2))
791         return;
792     }
793
794   /* Calculate results. */
795   for (i = 0; i < dsc->var_cnt; i++)
796     {
797       struct dsc_var *dv = &dsc->vars[i];
798       double W;
799       int j;
800
801       for (j = 0; j < DSC_N_STATS; j++)
802         dv->stats[j] = SYSMIS;
803
804       dv->valid = W = dsc->valid - dv->missing;
805
806       if (dv->moments != NULL)
807         moments_calculate (dv->moments, NULL,
808                            &dv->stats[DSC_MEAN], &dv->stats[DSC_VARIANCE],
809                            &dv->stats[DSC_SKEWNESS], &dv->stats[DSC_KURTOSIS]);
810       if (dsc->calc_stats & (1ul << DSC_SEMEAN)
811           && dv->stats[DSC_VARIANCE] != SYSMIS && W > 0.)
812         dv->stats[DSC_SEMEAN] = sqrt (dv->stats[DSC_VARIANCE]) / sqrt (W);
813       if (dsc->calc_stats & (1ul << DSC_STDDEV)
814           && dv->stats[DSC_VARIANCE] != SYSMIS)
815         dv->stats[DSC_STDDEV] = sqrt (dv->stats[DSC_VARIANCE]);
816       if (dsc->calc_stats & (1ul << DSC_SEKURT))
817         if (dv->stats[DSC_KURTOSIS] != SYSMIS)
818             dv->stats[DSC_SEKURT] = calc_sekurt (W);
819       if (dsc->calc_stats & (1ul << DSC_SESKEW)
820           && dv->stats[DSC_SKEWNESS] != SYSMIS)
821         dv->stats[DSC_SESKEW] = calc_seskew (W);
822       dv->stats[DSC_RANGE] = ((dv->min == DBL_MAX || dv->max == -DBL_MAX)
823                               ? SYSMIS : dv->max - dv->min);
824       dv->stats[DSC_MIN] = dv->min == DBL_MAX ? SYSMIS : dv->min;
825       dv->stats[DSC_MAX] = dv->max == -DBL_MAX ? SYSMIS : dv->max;
826       if (dsc->calc_stats & (1ul << DSC_SUM))
827         dv->stats[DSC_SUM] = W * dv->stats[DSC_MEAN];
828     }
829
830   /* Output results. */
831   display (dsc);
832 }
833
834 /* Returns true if any of the descriptives variables in DSC's
835    variable list have missing values in case C, false otherwise. */
836 static bool
837 listwise_missing (struct dsc_proc *dsc, const struct ccase *c)
838 {
839   size_t i;
840
841   for (i = 0; i < dsc->var_cnt; i++)
842     {
843       struct dsc_var *dv = &dsc->vars[i];
844       double x = case_num (c, dv->v);
845
846       if (var_is_num_missing (dv->v, x, dsc->exclude))
847         return true;
848     }
849   return false;
850 }
851 \f
852 /* Statistical display. */
853
854 static algo_compare_func descriptives_compare_dsc_vars;
855
856 /* Displays a table of descriptive statistics for DSC. */
857 static void
858 display (struct dsc_proc *dsc)
859 {
860   size_t i;
861   int nc;
862   struct tab_table *t;
863
864   nc = 1 + (dsc->format == DSC_SERIAL ? 2 : 1);
865   for (i = 0; i < DSC_N_STATS; i++)
866     if (dsc->show_stats & (1ul << i))
867       nc++;
868
869   if (dsc->sort_by_stat != DSC_NONE)
870     sort (dsc->vars, dsc->var_cnt, sizeof *dsc->vars,
871           descriptives_compare_dsc_vars, dsc);
872
873   t = tab_create (nc, dsc->var_cnt + 1);
874   tab_headers (t, 1, 0, 1, 0);
875   tab_box (t, TAL_1, TAL_1, -1, -1, 0, 0, nc - 1, dsc->var_cnt);
876   tab_box (t, -1, -1, -1, TAL_1, 1, 0, nc - 1, dsc->var_cnt);
877   tab_hline (t, TAL_2, 0, nc - 1, 1);
878   tab_vline (t, TAL_2, 1, 0, dsc->var_cnt);
879
880   nc = 0;
881   tab_text (t, nc++, 0, TAB_LEFT | TAT_TITLE, _("Variable"));
882   if (dsc->format == DSC_SERIAL)
883     {
884       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Valid N"));
885       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Missing N"));
886     }
887   else
888     tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, "N");
889
890   for (i = 0; i < DSC_N_STATS; i++)
891     if (dsc->show_stats & (1ul << i))
892       {
893         const char *title = gettext (dsc_info[i].name);
894         tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, title);
895       }
896
897   for (i = 0; i < dsc->var_cnt; i++)
898     {
899       struct dsc_var *dv = &dsc->vars[i];
900       size_t j;
901
902       nc = 0;
903       tab_text (t, nc++, i + 1, TAB_LEFT, var_get_name (dv->v));
904       tab_text_format (t, nc++, i + 1, 0, "%g", dv->valid);
905       if (dsc->format == DSC_SERIAL)
906         tab_text_format (t, nc++, i + 1, 0, "%g", dv->missing);
907
908       for (j = 0; j < DSC_N_STATS; j++)
909         if (dsc->show_stats & (1ul << j))
910           tab_double (t, nc++, i + 1, TAB_NONE, dv->stats[j], NULL);
911     }
912
913   tab_title (t, _("Valid cases = %g; cases with missing value(s) = %g."),
914              dsc->valid, dsc->missing_listwise);
915
916   tab_submit (t);
917 }
918
919 /* Compares `struct dsc_var's A and B according to the ordering
920    specified by CMD. */
921 static int
922 descriptives_compare_dsc_vars (const void *a_, const void *b_, const void *dsc_)
923 {
924   const struct dsc_var *a = a_;
925   const struct dsc_var *b = b_;
926   const struct dsc_proc *dsc = dsc_;
927
928   int result;
929
930   if (dsc->sort_by_stat == DSC_NAME)
931     result = strcasecmp (var_get_name (a->v), var_get_name (b->v));
932   else
933     {
934       double as = a->stats[dsc->sort_by_stat];
935       double bs = b->stats[dsc->sort_by_stat];
936
937       result = as < bs ? -1 : as > bs;
938     }
939
940   if (!dsc->sort_ascending)
941     result = -result;
942
943   return result;
944 }