Changed DFM from open-at-first-access to explicit-open. Before,
[pspp] / src / descript.q
1 /* PSPP - computes sample statistics.
2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3    Written by Ben Pfaff <blp@gnu.org>.
4
5    This program is free software; you can redistribute it and/or
6    modify it under the terms of the GNU General Public License as
7    published by the Free Software Foundation; either version 2 of the
8    License, or (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful, but
11    WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13    General Public License for more details.
14
15    You should have received a copy of the GNU General Public License
16    along with this program; if not, write to the Free Software
17    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
18    02111-1307, USA. */
19
20 /* FIXME: Many possible optimizations. */
21
22 #include <config.h>
23 #include <assert.h>
24 #include <limits.h>
25 #include <math.h>
26 #include <stdlib.h>
27 #include "algorithm.h"
28 #include "alloc.h"
29 #include "bitvector.h"
30 #include "command.h"
31 #include "lexer.h"
32 #include "error.h"
33 #include "magic.h"
34 #include "stats.h"
35 #include "som.h"
36 #include "tab.h"
37 #include "var.h"
38 #include "vfm.h"
39
40 /* (specification)
41    DESCRIPTIVES (dsc_):
42      *variables=custom;
43      +missing=miss:!variable/listwise,incl:!noinclude/include;
44      +format=labeled:!labels/nolabels,indexed:!noindex/index,lined:!line/serial;
45      +save=;
46      +options[op_]=1,2,3,4,5,6,7,8;
47      +statistics[st_]=all,1|mean,2|semean,5|stddev,6|variance,7|kurtosis,
48                       8|skewness,9|range,10|minimum,11|maximum,12|sum,
49                       13|default,14|seskewness,15|sekurtosis;
50      +sort=sortby:mean/semean/stddev/variance/kurtosis/skewness/range/
51            range/minimum/maximum/sum/name/seskewness/sekurtosis/!none, 
52            order:!a/d.
53 */
54 /* (declarations) */
55 /* (functions) */
56
57 /* DESCRIPTIVES private data. */
58
59 /* Describes properties of a distribution for the purpose of
60    calculating a Z-score. */
61 struct dsc_z_score
62   {
63     struct variable *s, *d;     /* Source, destination variable. */
64     double mean;                /* Distribution mean. */
65     double std_dev;             /* Distribution standard deviation. */
66   };
67
68 /* DESCRIPTIVES transformation (for calculating Z-scores). */
69 struct descriptives_trns
70   {
71     struct trns_header h;
72     int n;                      /* Number of Z-scores. */
73     struct dsc_z_score *z;      /* Array of Z-scores. */
74   };
75
76 /* These next three vars, see comment at top of display(). */
77 /* Number of cases missing listwise, even if option 5 not selected. */
78 static double d_glob_miss_list;
79
80 /* Number of total *cases* valid or missing, as a double.  Unless
81    option 5 is selected, d_glob_missing is 0. */
82 static double d_glob_valid, d_glob_missing;
83
84 /* Set when a weighting variable is missing or <=0. */
85 static int bad_weight;
86
87 /* Number of generic zvarnames we've generated in this execution. */
88 static int z_count;
89
90 /* Variables specified on command. */
91 static struct variable **v_variables;
92 static int n_variables;
93
94 /* Command specifications. */
95 static struct cmd_descriptives cmd;
96
97 /* Whether z-scores are computed. */
98 static int z_scores;
99
100 /* Statistic to sort by. */
101 static int sortby_stat;
102
103 /* Statistics to display. */
104 static unsigned long stats;
105
106 /* Easier access to long-named arrays. */
107 #define stat cmd.a_statistics
108 #define opt  cmd.a_options
109
110 /* Groups of statistics. */
111 #define BI          BIT_INDEX
112
113 #define dsc_default                                                     \
114         (BI (dsc_mean) | BI (dsc_stddev) | BI (dsc_min) | BI (dsc_max))
115      
116 #define dsc_all                                                 \
117         (BI (dsc_sum) | BI (dsc_min) | BI (dsc_max)             \
118          | BI (dsc_mean) | BI (dsc_semean) | BI (dsc_stddev)    \
119          | BI (dsc_variance) | BI (dsc_kurt) | BI (dsc_sekurt)  \
120          | BI (dsc_skew) | BI (dsc_seskew) | BI (dsc_range)     \
121          | BI (dsc_range))
122
123 /* Table of options. */
124 #define op_incl_miss    DSC_OP_1        /* Honored. */
125 #define op_no_varlabs   DSC_OP_2        /* Ignored. */
126 #define op_zscores      DSC_OP_3        /* Honored. */
127 #define op_index        DSC_OP_4        /* FIXME. */
128 #define op_excl_miss    DSC_OP_5        /* Honored. */
129 #define op_serial       DSC_OP_6        /* Honored. */
130 #define op_narrow       DSC_OP_7        /* Ignored. */
131 #define op_no_varnames  DSC_OP_8        /* Honored. */
132
133 /* Describes one statistic that can be calculated. */
134 /* FIXME: Currently sm,col_width are not used. */
135 struct dsc_info
136   {
137     int st_indx;                /* Index into st_a_statistics[]. */
138     int sb_indx;                /* Sort-by index. */
139     const char *s10;            /* Name, stuffed into 10 columns. */
140     const char *s8;             /* Name, stuffed into 8 columns. */
141     const char *sm;             /* Name, stuffed minimally. */
142     const char *s;              /* Full name. */
143     int max_degree;             /* Highest degree necessary to calculate this
144                                    statistic. */
145     int col_width;              /* Column width (not incl. spacing between columns) */
146   };
147
148 /* Table of statistics, indexed by dsc_*. */
149 static struct dsc_info dsc_info[dsc_n_stats] =
150 {
151   {DSC_ST_MEAN, DSC_MEAN, N_("Mean"), N_("Mean"), N_("Mean"), N_("mean"), 1, 10},
152   {DSC_ST_SEMEAN, DSC_SEMEAN, N_("S.E. Mean"), N_("S E Mean"), N_("SE"),
153    N_("standard error of mean"), 2, 10},
154   {DSC_ST_STDDEV, DSC_STDDEV, N_("Std Dev"), N_("Std Dev"), N_("SD"),
155    N_("standard deviation"), 2, 11},
156   {DSC_ST_VARIANCE, DSC_VARIANCE, N_("Variance"), N_("Variance"),
157    N_("Var"), N_("variance"), 2, 12},
158   {DSC_ST_KURTOSIS, DSC_KURTOSIS, N_("Kurtosis"), N_("Kurtosis"),
159    N_("Kurt"), N_("kurtosis"), 4, 9},
160   {DSC_ST_SEKURTOSIS, DSC_SEKURTOSIS, N_("S.E. Kurt"), N_("S E Kurt"), N_("SEKurt"),
161    N_("standard error of kurtosis"), 0, 9},
162   {DSC_ST_SKEWNESS, DSC_SKEWNESS, N_("Skewness"), N_("Skewness"), N_("Skew"),
163    N_("skewness"), 3, 9},
164   {DSC_ST_SESKEWNESS, DSC_SESKEWNESS, N_("S.E. Skew"), N_("S E Skew"), N_("SESkew"),
165    N_("standard error of skewness"), 0, 9},
166   {DSC_ST_RANGE, DSC_RANGE, N_("Range"), N_("Range"), N_("Rng"), N_("range"), 0, 10},
167   {DSC_ST_MINIMUM, DSC_MINIMUM, N_("Minimum"), N_("Minimum"), N_("Min"),
168    N_("minimum"), 0, 10},
169   {DSC_ST_MAXIMUM, DSC_MAXIMUM, N_("Maximum"), N_("Maximum"), N_("Max"),
170    N_("maximum"), 0, 10},
171   {DSC_ST_SUM, DSC_SUM, N_("Sum"), N_("Sum"), N_("Sum"), N_("sum"), 1, 13},
172 };
173
174 /* Z-score functions. */
175 static int generate_z_varname (struct variable * v);
176 static void dump_z_table (void);
177 static void run_z_pass (void);
178
179 /* Procedure execution functions. */
180 static int calc (struct ccase *, void *);
181 static void precalc (void *);
182 static void postcalc (void *);
183 static void display (void);
184 \f
185 /* Parser and outline. */
186
187 int
188 cmd_descriptives (void)
189 {
190   struct variable *v;
191   int i;
192
193   v_variables = NULL;
194   n_variables = 0;
195
196   if (!parse_descriptives (&cmd))
197     goto lossage;
198
199   if (n_variables == 0)
200     goto lossage;
201   for (i = 0; i < n_variables; i++)
202     {
203       v = v_variables[i];
204       v->p.dsc.dup = 0;
205       v->p.dsc.zname[0] = 0;
206     }
207
208   if (n_variables < 0)
209     {
210       msg (SE, _("No variables specified."));
211       goto lossage;
212     }
213
214   if (cmd.sbc_options && (cmd.sbc_save || cmd.sbc_format || cmd.sbc_missing))
215     {
216       msg (SE, _("OPTIONS may not be used with SAVE, FORMAT, or MISSING."));
217       goto lossage;
218     }
219   
220   if (!cmd.sbc_options)
221     {
222       if (cmd.incl == DSC_INCLUDE)
223         opt[op_incl_miss] = 1;
224       if (cmd.labeled == DSC_NOLABELS)
225         opt[op_no_varlabs] = 1;
226       if (cmd.sbc_save)
227         opt[op_zscores] = 1;
228       if (cmd.miss == DSC_LISTWISE)
229         opt[op_excl_miss] = 1;
230       if (cmd.lined == DSC_SERIAL)
231         opt[op_serial] = 1;
232     }
233
234   /* Construct z-score varnames, show translation table. */
235   if (opt[op_zscores])
236     {
237       z_count = 0;
238       for (i = 0; i < n_variables; i++)
239         {
240           v = v_variables[i];
241           if (v->p.dsc.dup++)
242             continue;
243
244           if (v->p.dsc.zname[0] == 0)
245             if (!generate_z_varname (v))
246               goto lossage;
247         }
248       dump_z_table ();
249       z_scores = 1;
250     }
251
252   /* Figure out statistics to calculate. */
253   stats = 0;
254   if (stat[DSC_ST_DEFAULT] || !cmd.sbc_statistics)
255     stats |= dsc_default;
256   if (stat[DSC_ST_ALL])
257     stats |= dsc_all;
258   for (i = 0; i < dsc_n_stats; i++)
259     if (stat[dsc_info[i].st_indx])
260       stats |= BIT_INDEX (i);
261   if (stats & dsc_kurt)
262     stats |= dsc_sekurt;
263   if (stats & dsc_skew)
264     stats |= dsc_seskew;
265
266   /* Check the sort order. */
267   sortby_stat = -1;
268   if (cmd.sortby == DSC_NONE)
269     sortby_stat = -2;
270   else if (cmd.sortby != DSC_NAME)
271     {
272       for (i = 0; i < n_variables; i++)
273         if (dsc_info[i].sb_indx == cmd.sortby)
274           {
275             sortby_stat = i;
276             if (!(stats & BIT_INDEX (i)))
277               {
278                 msg (SE, _("It's not possible to sort on `%s' without "
279                            "displaying `%s'."),
280                      gettext (dsc_info[i].s), gettext (dsc_info[i].s));
281                 goto lossage;
282               }
283           }
284       assert (sortby_stat != -1);
285     }
286
287   /* Data pass! */
288   bad_weight = 0;
289   procedure_with_splits (precalc, calc, postcalc, NULL);
290
291   if (bad_weight)
292     msg (SW, _("At least one case in the data file had a weight value "
293          "that was system-missing, zero, or negative.  These case(s) "
294          "were ignored."));
295
296   /* Z-scoring! */
297   if (z_scores)
298     run_z_pass ();
299
300   if (v_variables)
301     free (v_variables);
302   return CMD_SUCCESS;
303
304  lossage:
305   if (v_variables)
306     free (v_variables);
307   return CMD_FAILURE;
308 }
309
310 /* Parses the VARIABLES subcommand. */
311 static int
312 dsc_custom_variables (struct cmd_descriptives *cmd UNUSED)
313 {
314   if (!lex_match_id ("VARIABLES")
315       && (token != T_ID || dict_lookup_var (default_dict, tokid) == NULL)
316       && token != T_ALL)
317     return 2;
318   lex_match ('=');
319
320   while (token == T_ID || token == T_ALL)
321     {
322       int i, n;
323
324       n = n_variables;
325       if (!parse_variables (default_dict, &v_variables, &n_variables,
326                             PV_DUPLICATE | PV_SINGLE | PV_APPEND | PV_NUMERIC
327                             | PV_NO_SCRATCH))
328         return 0;
329       if (lex_match ('('))
330         {
331           if (n_variables - n > 1)
332             {
333               msg (SE, _("Names for z-score variables must be given for "
334                          "individual variables, not for groups of "
335                          "variables."));
336               return 0;
337             }
338           assert (n_variables - n <= 0);
339           if (token != T_ID)
340             {
341               msg (SE, _("Name for z-score variable expected."));
342               return 0;
343             }
344           if (dict_lookup_var (default_dict, tokid))
345             {
346               msg (SE, _("Z-score variable name `%s' is a "
347                          "duplicate variable name with a current variable."),
348                    tokid);
349               return 0;
350             }
351           for (i = 0; i < n_variables; i++)
352             if (v_variables[i]->p.dsc.zname[0]
353                 && !strcmp (v_variables[i]->p.dsc.zname, tokid))
354               {
355                 msg (SE, _("Z-score variable name `%s' is "
356                            "used multiple times."), tokid);
357                 return 0;
358               }
359           strcpy (v_variables[n_variables - 1]->p.dsc.zname, tokid);
360           lex_get ();
361           if (token != ')')
362             {
363               msg (SE, _("`)' expected after z-score variable name."));
364               return 0;
365             }
366
367           z_scores = 1;
368         }
369       lex_match (',');
370     }
371   return 1;
372 }
373 \f
374 /* Z scores. */
375
376 /* Returns 0 if NAME is a duplicate of any existing variable name or
377    of any previously-declared z-var name; otherwise returns 1. */
378 static int
379 try_name (char *name)
380 {
381   int i;
382
383   if (dict_lookup_var (default_dict, name) != NULL)
384     return 0;
385   for (i = 0; i < n_variables; i++)
386     {
387       struct variable *v = v_variables[i];
388       if (!strcmp (v->p.dsc.zname, name))
389         return 0;
390     }
391   return 1;
392 }
393
394 static int
395 generate_z_varname (struct variable * v)
396 {
397   char zname[10];
398
399   strcpy (&zname[1], v->name);
400   zname[0] = 'Z';
401   zname[8] = '\0';
402   if (try_name (zname))
403     {
404       strcpy (v->p.dsc.zname, zname);
405       return 1;
406     }
407
408   for (;;)
409     {
410       /* Generate variable name. */
411       z_count++;
412
413       if (z_count <= 99)
414         sprintf (zname, "ZSC%03d", z_count);
415       else if (z_count <= 108)
416         sprintf (zname, "STDZ%02d", z_count - 99);
417       else if (z_count <= 117)
418         sprintf (zname, "ZZZZ%02d", z_count - 108);
419       else if (z_count <= 126)
420         sprintf (zname, "ZQZQ%02d", z_count - 117);
421       else
422         {
423           msg (SE, _("Ran out of generic names for Z-score variables.  "
424                      "There are only 126 generic names: ZSC001-ZSC0999, "
425                      "STDZ01-STDZ09, ZZZZ01-ZZZZ09, ZQZQ01-ZQZQ09."));
426           return 0;
427         }
428       
429       if (try_name (zname))
430         {
431           strcpy (v->p.dsc.zname, zname);
432           return 1;
433         }
434     }
435 }
436
437 static void
438 dump_z_table (void)
439 {
440   int count;
441   struct tab_table *t;
442   
443   {
444     int i;
445     
446     for (count = i = 0; i < n_variables; i++)
447       if (v_variables[i]->p.dsc.zname)
448         count++;
449   }
450   
451   t = tab_create (2, count + 1, 0);
452   tab_title (t, 0, _("Mapping of variables to corresponding Z-scores."));
453   tab_columns (t, SOM_COL_DOWN, 1);
454   tab_headers (t, 0, 0, 1, 0);
455   tab_box (t, TAL_1, TAL_1, TAL_0, TAL_1, 0, 0, 1, count);
456   tab_hline (t, TAL_2, 0, 1, 1);
457   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
458   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Target"));
459   tab_dim (t, tab_natural_dimensions);
460
461   {
462     int i, y;
463     
464     for (i = 0, y = 1; i < n_variables; i++)
465       if (v_variables[i]->p.dsc.zname)
466         {
467           tab_text (t, 0, y, TAB_LEFT, v_variables[i]->name);
468           tab_text (t, 1, y++, TAB_LEFT, v_variables[i]->p.dsc.zname);
469         }
470   }
471   
472   tab_submit (t);
473 }
474
475 /* Transformation function to calculate Z-scores. */
476 static int
477 descriptives_trns_proc (struct trns_header * trns, struct ccase * c,
478                         int case_num UNUSED)
479 {
480   struct descriptives_trns *t = (struct descriptives_trns *) trns;
481   struct dsc_z_score *z;
482   int i;
483
484   for (i = t->n, z = t->z; i--; z++)
485     {
486       double score = c->data[z->s->fv].f;
487
488       if (z->mean == SYSMIS || score == SYSMIS)
489         c->data[z->d->fv].f = SYSMIS;
490       else
491         c->data[z->d->fv].f = (score - z->mean) / z->std_dev;
492     }
493   return -1;
494 }
495
496 /* Frees a descriptives_trns struct. */
497 static void
498 descriptives_trns_free (struct trns_header * trns)
499 {
500   struct descriptives_trns *t = (struct descriptives_trns *) trns;
501
502   free (t->z);
503 }
504
505 /* The name is a misnomer: actually this function sets up a
506    transformation by which scores can be converted into Z-scores. */
507 static void
508 run_z_pass (void)
509 {
510   struct descriptives_trns *t;
511   int count, i;
512
513   for (i = 0; i < n_variables; i++)
514     v_variables[i]->p.dsc.dup = 0;
515   for (count = i = 0; i < n_variables; i++)
516     {
517       if (v_variables[i]->p.dsc.dup++)
518         continue;
519       if (v_variables[i]->p.dsc.zname)
520         count++;
521     }
522
523   t = xmalloc (sizeof *t);
524   t->h.proc = descriptives_trns_proc;
525   t->h.free = descriptives_trns_free;
526   t->n = count;
527   t->z = xmalloc (count * sizeof *t->z);
528
529   for (i = 0; i < n_variables; i++)
530     v_variables[i]->p.dsc.dup = 0;
531   for (count = i = 0; i < n_variables; i++)
532     {
533       struct variable *v = v_variables[i];
534       if (v->p.dsc.dup++ == 0 && v->p.dsc.zname[0])
535         {
536           char *cp;
537           struct variable *d;
538
539           t->z[count].s = v;
540           t->z[count].d = d = dict_create_var_assert (default_dict,
541                                                       v->p.dsc.zname, 0);
542           d->init = 0;
543           if (v->label)
544             {
545               d->label = xmalloc (strlen (v->label) + 12);
546               cp = stpcpy (d->label, _("Z-score of "));
547               strcpy (cp, v->label);
548             }
549           else
550             {
551               d->label = xmalloc (strlen (v->name) + 12);
552               cp = stpcpy (d->label, _("Z-score of "));
553               strcpy (cp, v->name);
554             }
555           t->z[count].mean = v->p.dsc.stats[dsc_mean];
556           t->z[count].std_dev = v->p.dsc.stats[dsc_stddev];
557           if (t->z[count].std_dev == SYSMIS || t->z[count].std_dev == 0.0)
558             t->z[count].mean = SYSMIS;
559           count++;
560         }
561     }
562
563   add_transformation ((struct trns_header *) t);
564 }
565 \f
566 /* Statistical calculation. */
567
568 static void
569 precalc (void *aux UNUSED)
570 {
571   int i;
572
573   for (i = 0; i < n_variables; i++)
574     v_variables[i]->p.dsc.dup = -2;
575   for (i = 0; i < n_variables; i++)
576     {
577       struct descriptives_proc *dsc = &v_variables[i]->p.dsc;
578
579       /* Don't need to initialize more than once. */
580       if (dsc->dup == -1)
581         continue;
582       dsc->dup = -1;
583
584       dsc->valid = dsc->miss = 0.0;
585       dsc->X_bar = dsc->M2 = dsc->M3 = dsc->M4 = 0.0;
586       dsc->min = DBL_MAX;
587       dsc->max = -DBL_MAX;
588     }
589   d_glob_valid = d_glob_missing = 0.0;
590 }
591
592 static int
593 calc (struct ccase * c, void *aux UNUSED)
594 {
595   int i;
596
597   /* Unique case identifier. */
598   static int case_id;
599
600   /* Get the weight for this case. */
601   double weight = dict_get_case_weight (default_dict, c);
602
603   if (weight <= 0.0)
604     {
605       weight = 0.0;
606       bad_weight = 1;
607     }
608   case_id++;
609
610   /* Handle missing values. */
611   for (i = 0; i < n_variables; i++)
612     {
613       struct variable *v = v_variables[i];
614       double X = c->data[v->fv].f;
615
616       if (X == SYSMIS || (!opt[op_incl_miss] && is_num_user_missing (X, v)))
617         {
618           if (opt[op_excl_miss])
619             {
620               d_glob_missing += weight;
621               return 1;
622             }
623           else
624             {
625               d_glob_miss_list += weight;
626               goto iterate;
627             }
628         }
629     }
630   d_glob_valid += weight;
631
632 iterate:
633   for (i = 0; i < n_variables; i++)
634     {
635       struct descriptives_proc *inf = &v_variables[i]->p.dsc;
636
637       double X, v;
638       double W_old, W_new;
639       double v2, v3, v4;
640       double w2, w3, w4;
641
642       if (inf->dup == case_id)
643         continue;
644       inf->dup = case_id;
645
646       X = c->data[v_variables[i]->fv].f;
647       if (X == SYSMIS
648           || (!opt[op_incl_miss] && is_num_user_missing (X, v_variables[i])))
649         {
650           inf->miss += weight;
651           continue;
652         }
653
654       /* These formulas taken from _SPSS Statistical Algorithms_.  The
655          names W_old, and W_new are used for W_j-1 and W_j,
656          respectively, and other variables simply have the subscripts
657          trimmed off, except for X_bar.
658
659          I am happy that mathematical formulas may not be
660          copyrighted. */
661       W_old = inf->valid;
662       W_new = inf->valid += weight;
663       v = (weight / W_new) * (X - inf->X_bar);
664       v2 = v * v;
665       v3 = v2 * v;
666       v4 = v3 * v;
667       w2 = weight * weight;
668       w3 = w2 * weight;
669       w4 = w3 * weight;
670       inf->M4 += (-4.0 * v * inf->M3 + 6.0 * v2 * inf->M2
671                + (W_new * W_new - 3 * weight * W_old / w3) * v4 * W_old * W_new);
672       inf->M3 += (-3.0 * v * inf->M2 + W_new * W_old / w2
673                   * (W_new - 2 * weight) * v3);
674       inf->M2 += W_new * W_old / weight * v2;
675       inf->X_bar += v;
676       if (X < inf->min)
677         inf->min = X;
678       if (X > inf->max)
679         inf->max = X;
680     }
681   return 1;
682 }
683
684 static void
685 postcalc (void *aux UNUSED)
686 {
687   int i;
688
689   if (opt[op_excl_miss])
690     d_glob_miss_list = d_glob_missing;
691
692   for (i = 0; i < n_variables; i++)
693     {
694       struct descriptives_proc *dsc = &v_variables[i]->p.dsc;
695       double W;
696
697       /* Don't duplicate our efforts. */
698       if (dsc->dup == -2)
699         continue;
700       dsc->dup = -2;
701
702       W = dsc->valid;
703
704       dsc->stats[dsc_mean] = dsc->X_bar;
705       dsc->stats[dsc_variance] = dsc->M2 / (W - 1);
706       dsc->stats[dsc_stddev] = sqrt (dsc->stats[dsc_variance]);
707       dsc->stats[dsc_semean] = dsc->stats[dsc_stddev] / sqrt (W);
708       dsc->stats[dsc_min] = dsc->min == DBL_MAX ? SYSMIS : dsc->min;
709       dsc->stats[dsc_max] = dsc->max == -DBL_MAX ? SYSMIS : dsc->max;
710       dsc->stats[dsc_range] = ((dsc->min == DBL_MAX || dsc->max == -DBL_MAX)
711                                ? SYSMIS : dsc->max - dsc->min);
712       dsc->stats[dsc_sum] = W * dsc->X_bar;
713       if (W > 2.0 && dsc->stats[dsc_variance] >= 1e-20)
714         {
715           double S = dsc->stats[dsc_stddev];
716           dsc->stats[dsc_skew] = (W * dsc->M3 / ((W - 1.0) * (W - 2.0) * S * S * S));
717           dsc->stats[dsc_seskew] =
718             sqrt (6.0 * W * (W - 1.0) / ((W - 2.0) * (W + 1.0) * (W + 3.0)));
719         }
720       else
721         {
722           dsc->stats[dsc_skew] = dsc->stats[dsc_seskew] = SYSMIS;
723         }
724       if (W > 3.0 && dsc->stats[dsc_variance] >= 1e-20)
725         {
726           double S2 = dsc->stats[dsc_variance];
727           double SE_g1 = dsc->stats[dsc_seskew];
728
729           dsc->stats[dsc_kurt] =
730             (W * (W + 1.0) * dsc->M4 - 3.0 * dsc->M2 * dsc->M2 * (W - 1.0))
731             / ((W - 1.0) * (W - 2.0) * (W - 3.0) * S2 * S2);
732
733           /* Note that in _SPSS Statistical Algorithms_, the square
734              root symbol is missing from this formula. */
735           dsc->stats[dsc_sekurt] =
736             sqrt ((4.0 * (W * W - 1.0) * SE_g1 * SE_g1) / ((W - 3.0) * (W + 5.0)));
737         }
738       else
739         {
740           dsc->stats[dsc_kurt] = dsc->stats[dsc_sekurt] = SYSMIS;
741         }
742     }
743
744   display ();
745 }
746 \f
747 /* Statistical display. */
748
749 static algo_compare_func descriptives_compare_variables;
750
751 static void
752 display (void)
753 {
754   int i, j;
755
756   int nc, n_stats;
757   struct tab_table *t;
758
759   /* If op_excl_miss is on, d_glob_valid and (potentially)
760      d_glob_missing are nonzero, and d_glob_missing equals
761      d_glob_miss_list.
762
763      If op_excl_miss is off, d_glob_valid is nonzero.  d_glob_missing
764      is zero.  d_glob_miss_list is (potentially) nonzero.  */
765
766   if (sortby_stat != -2)
767     sort (v_variables, n_variables, sizeof *v_variables,
768            descriptives_compare_variables, &cmd);
769
770   for (nc = i = 0; i < dsc_n_stats; i++)
771     if (stats & BIT_INDEX (i))
772       nc++;
773   n_stats = nc;
774   if (!opt[op_no_varnames])
775     nc++;
776   nc += opt[op_serial] ? 2 : 1;
777
778   t = tab_create (nc, n_variables + 1, 0);
779   tab_headers (t, 1, 0, 1, 0);
780   tab_box (t, TAL_1, TAL_1, -1, -1, 0, 0, nc - 1, n_variables);
781   tab_box (t, -1, -1, -1, TAL_1, 1, 0, nc - 1, n_variables);
782   tab_hline (t, TAL_2, 0, nc - 1, 1);
783   tab_vline (t, TAL_2, 1, 0, n_variables);
784   tab_dim (t, tab_natural_dimensions);
785
786   nc = 0;
787   if (!opt[op_no_varnames])
788     {
789       tab_text (t, nc++, 0, TAB_LEFT | TAT_TITLE, _("Variable"));
790     }
791   if (opt[op_serial])
792     {
793       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Valid N"));
794       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Missing N"));
795     } else {
796       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, "N");
797     }
798
799   for (i = 0; i < dsc_n_stats; i++)
800     if (stats & BIT_INDEX (i))
801       {
802         const char *title = gettext (dsc_info[i].s8);
803         tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, title);
804       }
805
806   for (i = 0; i < n_variables; i++)
807     {
808       struct variable *v = v_variables[i];
809
810       nc = 0;
811       if (!opt[op_no_varnames])
812         tab_text (t, nc++, i + 1, TAB_LEFT, v->name);
813       tab_text (t, nc++, i + 1, TAT_PRINTF, "%g", v->p.dsc.valid);
814       if (opt[op_serial])
815         tab_text (t, nc++, i + 1, TAT_PRINTF, "%g", v->p.dsc.miss);
816       for (j = 0; j < dsc_n_stats; j++)
817         if (stats & BIT_INDEX (j))
818           tab_float (t, nc++, i + 1, TAB_NONE, v->p.dsc.stats[j], 10, 3);
819     }
820
821   tab_title (t, 1, _("Valid cases = %g; cases with missing value(s) = %g."),
822              d_glob_valid, d_glob_miss_list);
823
824   tab_submit (t);
825 }
826
827 /* Compares variables A and B according to the ordering specified
828    by CMD. */
829 static int
830 descriptives_compare_variables (const void *a_, const void *b_, void *cmd_)
831 {
832   struct variable *const *ap = a_;
833   struct variable *const *bp = b_;
834   struct variable *a = *ap;
835   struct variable *b = *bp;
836   struct cmd_descriptives *cmd = cmd_;
837
838   int result;
839
840   if (cmd->sortby == DSC_NAME)
841     result = strcmp (a->name, b->name);
842   else 
843     {
844       double as = a->p.dsc.stats[sortby_stat];
845       double bs = b->p.dsc.stats[sortby_stat];
846
847       result = as < bs ? -1 : as > bs;
848     }
849
850   if (cmd->order == DSC_D)
851     result = -result;
852
853   return result;
854 }
855
856 /*
857    Local variables:
858    mode: c
859    End:
860 */