src/language/stats/descriptives.c

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