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