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
  24 #include <limits.h>
  25 #include <math.h>
  26 #include <stdlib.h>
  27
  28 #include <data/case.h>
  29 #include <data/casefile.h>
  30 #include <data/dictionary.h>
  31 #include <data/procedure.h>
  32 #include <data/transformations.h>
  33 #include <data/variable.h>
  34 #include <language/command.h>
  35 #include <language/dictionary/split-file.h>
  36 #include <language/lexer/lexer.h>
  37 #include <language/lexer/variable-parser.h>
  38 #include <libpspp/alloc.h>
  39 #include <libpspp/array.h>
  40 #include <libpspp/compiler.h>
  41 #include <libpspp/magic.h>
  42 #include <libpspp/message.h>
  43 #include <libpspp/assertion.h>
  44 #include <math/moments.h>
  45 #include <output/manager.h>
  46 #include <output/table.h>
  47
  48 #include "gettext.h"
  49 #define _(msgid) gettext (msgid)
  50 #define N_(msgid) msgid
  51
  52 /* DESCRIPTIVES private data. */
  53
  54 struct dsc_proc;
  55
  56 /* Handling of missing values. */
  57 enum dsc_missing_type
  58   {
  59     DSC_VARIABLE,       /* Handle missing values on a per-variable basis. */
  60     DSC_LISTWISE        /* Discard entire case if any variable is missing. */
  61   };
  62
  63 /* Describes properties of a distribution for the purpose of
  64    calculating a Z-score. */
  65 struct dsc_z_score
  66   {
  67     int src_idx;                /* Source index into case data. */
  68     int dst_idx;                /* Destination index into case data. */
  69     double mean;                /* Distribution mean. */
  70     double std_dev;             /* Distribution standard deviation. */
  71     struct variable *v;         /* Variable on which z-score is based. */
  72   };
  73
  74 /* DESCRIPTIVES transformation (for calculating Z-scores). */
  75 struct dsc_trns
  76   {
  77     struct dsc_z_score *z_scores; /* Array of Z-scores. */
  78     int z_score_cnt;            /* Number of Z-scores. */
  79     struct variable **vars;     /* Variables for listwise missing checks. */
  80     size_t var_cnt;             /* Number of variables. */
  81     enum dsc_missing_type missing_type; /* Treatment of missing values. */
  82     int include_user_missing;   /* Nonzero to include user-missing values. */
  83   };
  84
  85 /* Statistics.  Used as bit indexes, so must be 32 or fewer. */
  86 enum dsc_statistic
  87   {
  88     DSC_MEAN = 0, DSC_SEMEAN, DSC_STDDEV, DSC_VARIANCE, DSC_KURTOSIS,
  89     DSC_SEKURT, DSC_SKEWNESS, DSC_SESKEW, DSC_RANGE, DSC_MIN,
  90     DSC_MAX, DSC_SUM, DSC_N_STATS,
  91
  92     /* Only valid as sort criteria. */
  93     DSC_NAME = -2,              /* Sort by name. */
  94     DSC_NONE = -1               /* Unsorted. */
  95   };
  96
  97 /* Describes one statistic. */
  98 struct dsc_statistic_info
  99   {
 100     const char *identifier;     /* Identifier. */
 101     const char *name;           /* Full name. */
 102     enum moment moment;         /* Highest moment needed to calculate. */
 103   };
 104
 105 /* Table of statistics, indexed by DSC_*. */
 106 static const struct dsc_statistic_info dsc_info[DSC_N_STATS] =
 107   {
 108     {"MEAN", N_("Mean"), MOMENT_MEAN},
 109     {"SEMEAN", N_("S E Mean"), MOMENT_VARIANCE},
 110     {"STDDEV", N_("Std Dev"), MOMENT_VARIANCE},
 111     {"VARIANCE", N_("Variance"), MOMENT_VARIANCE},
 112     {"KURTOSIS", N_("Kurtosis"), MOMENT_KURTOSIS},
 113     {"SEKURTOSIS", N_("S E Kurt"), MOMENT_NONE},
 114     {"SKEWNESS", N_("Skewness"), MOMENT_SKEWNESS},
 115     {"SESKEWNESS", N_("S E Skew"), MOMENT_NONE},
 116     {"RANGE", N_("Range"), MOMENT_NONE},
 117     {"MINIMUM", N_("Minimum"), MOMENT_NONE},
 118     {"MAXIMUM", N_("Maximum"), MOMENT_NONE},
 119     {"SUM", N_("Sum"), MOMENT_MEAN},
 120   };
 121
 122 /* Statistics calculated by default if none are explicitly
 123    requested. */
 124 #define DEFAULT_STATS                                                   \
 125         ((1ul << DSC_MEAN) | (1ul << DSC_STDDEV) | (1ul << DSC_MIN)     \
 126          | (1ul << DSC_MAX))
 127
 128 /* A variable specified on DESCRIPTIVES. */
 129 struct dsc_var
 130   {
 131     struct variable *v;         /* Variable to calculate on. */
 132     char z_name[LONG_NAME_LEN + 1]; /* Name for z-score variable. */
 133     double valid, missing;      /* Valid, missing counts. */
 134     struct moments *moments;    /* Moments. */
 135     double min, max;            /* Maximum and mimimum values. */
 136     double stats[DSC_N_STATS];  /* All the stats' values. */
 137   };
 138
 139 /* Output format. */
 140 enum dsc_format
 141   {
 142     DSC_LINE,           /* Abbreviated format. */
 143     DSC_SERIAL          /* Long format. */
 144   };
 145
 146 /* A DESCRIPTIVES procedure. */
 147 struct dsc_proc
 148   {
 149     /* Per-variable info. */
 150     struct dsc_var *vars;       /* Variables. */
 151     size_t var_cnt;             /* Number of variables. */
 152
 153     /* User options. */
 154     enum dsc_missing_type missing_type; /* Treatment of missing values. */
 155     int include_user_missing;   /* Nonzero to include user-missing values. */
 156     int show_var_labels;        /* Nonzero to show variable labels. */
 157     int show_index;             /* Nonzero to show variable index. */
 158     enum dsc_format format;     /* Output format. */
 159
 160     /* Accumulated results. */
 161     double missing_listwise;    /* Sum of weights of cases missing listwise. */
 162     double valid;               /* Sum of weights of valid cases. */
 163     bool bad_warn;               /* Warn if bad weight found. */
 164     enum dsc_statistic sort_by_stat; /* Statistic to sort by; -1: name. */
 165     int sort_ascending;         /* !0: ascending order; 0: descending. */
 166     unsigned long show_stats;   /* Statistics to display. */
 167     unsigned long calc_stats;   /* Statistics to calculate. */
 168     enum moment max_moment;     /* Highest moment needed for stats. */
 169   };
 170
 171 /* Parsing. */
 172 static enum dsc_statistic match_statistic (struct lexer *);
 173 static void free_dsc_proc (struct dsc_proc *);
 174
 175 /* Z-score functions. */
 176 static bool try_name (const struct dictionary *dict,
 177                       struct dsc_proc *dsc, const char *name);
 178 static bool generate_z_varname (const struct dictionary *dict,
 179                                 struct dsc_proc *dsc, char *z_name,
 180                                 const char *name, size_t *z_cnt);
 181 static void dump_z_table (struct dsc_proc *);
 182 static void setup_z_trns (struct dsc_proc *, struct dataset *);
 183
 184 /* Procedure execution functions. */
 185 static bool calc_descriptives (const struct ccase *first,
 186                                const struct casefile *, void *dsc_,
 187                                const struct dataset *);
 188 static void display (struct dsc_proc *dsc);
 189 \f
 190 /* Parser and outline. */
 191
 192 /* Handles DESCRIPTIVES. */
 193 int
 194 cmd_descriptives (struct lexer *lexer, struct dataset *ds)
 195 {
 196   struct dictionary *dict = dataset_dict (ds);
 197   struct dsc_proc *dsc;
 198   struct variable **vars = NULL;
 199   size_t var_cnt = 0;
 200   int save_z_scores = 0;
 201   size_t z_cnt = 0;
 202   size_t i;
 203   bool ok;
 204
 205   /* Create and initialize dsc. */
 206   dsc = xmalloc (sizeof *dsc);
 207   dsc->vars = NULL;
 208   dsc->var_cnt = 0;
 209   dsc->missing_type = DSC_VARIABLE;
 210   dsc->include_user_missing = 0;
 211   dsc->show_var_labels = 1;
 212   dsc->show_index = 0;
 213   dsc->format = DSC_LINE;
 214   dsc->missing_listwise = 0.;
 215   dsc->valid = 0.;
 216   dsc->bad_warn = 1;
 217   dsc->sort_by_stat = DSC_NONE;
 218   dsc->sort_ascending = 1;
 219   dsc->show_stats = dsc->calc_stats = DEFAULT_STATS;
 220
 221   /* Parse DESCRIPTIVES. */
 222   while (lex_token (lexer) != '.')
 223     {
 224       if (lex_match_id (lexer, "MISSING"))
 225         {
 226           lex_match (lexer, '=');
 227           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 228             {
 229               if (lex_match_id (lexer, "VARIABLE"))
 230                 dsc->missing_type = DSC_VARIABLE;
 231               else if (lex_match_id (lexer, "LISTWISE"))
 232                 dsc->missing_type = DSC_LISTWISE;
 233               else if (lex_match_id (lexer, "INCLUDE"))
 234                 dsc->include_user_missing = 1;
 235               else
 236                 {
 237                   lex_error (lexer, NULL);
 238                   goto error;
 239                 }
 240               lex_match (lexer, ',');
 241             }
 242         }
 243       else if (lex_match_id (lexer, "SAVE"))
 244         save_z_scores = 1;
 245       else if (lex_match_id (lexer, "FORMAT"))
 246         {
 247           lex_match (lexer, '=');
 248           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 249             {
 250               if (lex_match_id (lexer, "LABELS"))
 251                 dsc->show_var_labels = 1;
 252               else if (lex_match_id (lexer, "NOLABELS"))
 253                 dsc->show_var_labels = 0;
 254               else if (lex_match_id (lexer, "INDEX"))
 255                 dsc->show_index = 1;
 256               else if (lex_match_id (lexer, "NOINDEX"))
 257                 dsc->show_index = 0;
 258               else if (lex_match_id (lexer, "LINE"))
 259                 dsc->format = DSC_LINE;
 260               else if (lex_match_id (lexer, "SERIAL"))
 261                 dsc->format = DSC_SERIAL;
 262               else
 263                 {
 264                   lex_error (lexer, NULL);
 265                   goto error;
 266                 }
 267               lex_match (lexer, ',');
 268             }
 269         }
 270       else if (lex_match_id (lexer, "STATISTICS"))
 271         {
 272           lex_match (lexer, '=');
 273           dsc->show_stats = 0;
 274           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 275             {
 276               if (lex_match (lexer, T_ALL))
 277                 dsc->show_stats |= (1ul << DSC_N_STATS) - 1;
 278               else if (lex_match_id (lexer, "DEFAULT"))
 279                 dsc->show_stats |= DEFAULT_STATS;
 280               else
 281                 dsc->show_stats |= 1ul << (match_statistic (lexer));
 282               lex_match (lexer, ',');
 283             }
 284           if (dsc->show_stats == 0)
 285             dsc->show_stats = DEFAULT_STATS;
 286         }
 287       else if (lex_match_id (lexer, "SORT"))
 288         {
 289           lex_match (lexer, '=');
 290           if (lex_match_id (lexer, "NAME"))
 291             dsc->sort_by_stat = DSC_NAME;
 292           else
 293             {
 294               dsc->sort_by_stat = match_statistic (lexer);
 295               if (dsc->sort_by_stat == DSC_NONE )
 296                 dsc->sort_by_stat = DSC_MEAN;
 297             }
 298           if (lex_match (lexer, '('))
 299             {
 300               if (lex_match_id (lexer, "A"))
 301                 dsc->sort_ascending = 1;
 302               else if (lex_match_id (lexer, "D"))
 303                 dsc->sort_ascending = 0;
 304               else
 305                 lex_error (lexer, NULL);
 306               lex_force_match (lexer, ')');
 307             }
 308         }
 309       else if (var_cnt == 0)
 310         {
 311           if (lex_look_ahead (lexer) == '=')
 312             {
 313               lex_match_id (lexer, "VARIABLES");
 314               lex_match (lexer, '=');
 315             }
 316
 317           while (lex_token (lexer) != '.' && lex_token (lexer) != '/')
 318             {
 319               int i;
 320
 321               if (!parse_variables (lexer, dataset_dict (ds), &vars, &var_cnt,
 322                                     PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
 323                 goto error;
 324
 325               dsc->vars = xnrealloc (dsc->vars, var_cnt, sizeof *dsc->vars);
 326               for (i = dsc->var_cnt; i < var_cnt; i++)
 327                 {
 328                   struct dsc_var *dv = &dsc->vars[i];
 329                   dv->v = vars[i];
 330                   dv->z_name[0] = '\0';
 331                   dv->moments = NULL;
 332                 }
 333               dsc->var_cnt = var_cnt;
 334
 335               if (lex_match (lexer, '('))
 336                 {
 337                   if (lex_token (lexer) != T_ID)
 338                     {
 339                       lex_error (lexer, NULL);
 340                       goto error;
 341                     }
 342                   if (try_name (dict, dsc, lex_tokid (lexer)))
 343                     {
 344                       strcpy (dsc->vars[dsc->var_cnt - 1].z_name, lex_tokid (lexer));
 345                       z_cnt++;
 346                     }
 347                   else
 348                     msg (SE, _("Z-score variable name %s would be"
 349                                " a duplicate variable name."), lex_tokid (lexer));
 350                   lex_get (lexer);
 351                   if (!lex_force_match (lexer, ')'))
 352                     goto error;
 353                 }
 354             }
 355         }
 356       else
 357         {
 358           lex_error (lexer, NULL);
 359           goto error;
 360         }
 361
 362       lex_match (lexer, '/');
 363     }
 364   if (var_cnt == 0)
 365     {
 366       msg (SE, _("No variables specified."));
 367       goto error;
 368     }
 369
 370   /* Construct z-score varnames, show translation table. */
 371   if (z_cnt || save_z_scores)
 372     {
 373       if (save_z_scores)
 374         {
 375           size_t gen_cnt = 0;
 376
 377           for (i = 0; i < dsc->var_cnt; i++)
 378             if (dsc->vars[i].z_name[0] == 0)
 379               {
 380                 if (!generate_z_varname (dict, dsc, dsc->vars[i].z_name,
 381                                          var_get_name (dsc->vars[i].v),
 382                                          &gen_cnt))
 383                   goto error;
 384                 z_cnt++;
 385               }
 386         }
 387       dump_z_table (dsc);
 388     }
 389
 390   /* Figure out statistics to display. */
 391   if (dsc->show_stats & (1ul << DSC_SKEWNESS))
 392     dsc->show_stats |= 1ul << DSC_SESKEW;
 393   if (dsc->show_stats & (1ul << DSC_KURTOSIS))
 394     dsc->show_stats |= 1ul << DSC_SEKURT;
 395
 396   /* Figure out which statistics to calculate. */
 397   dsc->calc_stats = dsc->show_stats;
 398   if (z_cnt > 0)
 399     dsc->calc_stats |= (1ul << DSC_MEAN) | (1ul << DSC_STDDEV);
 400   if (dsc->sort_by_stat >= 0)
 401     dsc->calc_stats |= 1ul << dsc->sort_by_stat;
 402   if (dsc->show_stats & (1ul << DSC_SESKEW))
 403     dsc->calc_stats |= 1ul << DSC_SKEWNESS;
 404   if (dsc->show_stats & (1ul << DSC_SEKURT))
 405     dsc->calc_stats |= 1ul << DSC_KURTOSIS;
 406
 407   /* Figure out maximum moment needed and allocate moments for
 408      the variables. */
 409   dsc->max_moment = MOMENT_NONE;
 410   for (i = 0; i < DSC_N_STATS; i++)
 411     if (dsc->calc_stats & (1ul << i) && dsc_info[i].moment > dsc->max_moment)
 412       dsc->max_moment = dsc_info[i].moment;
 413   if (dsc->max_moment != MOMENT_NONE)
 414     for (i = 0; i < dsc->var_cnt; i++)
 415       dsc->vars[i].moments = moments_create (dsc->max_moment);
 416
 417   /* Data pass. */
 418   ok = multipass_procedure_with_splits (ds, calc_descriptives, dsc);
 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, size_t *z_cnt)
 497 {
 498   char name[LONG_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);
 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   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)->fv);
 600           if ( score == SYSMIS
 601                || (!t->include_user_missing
 602                    && var_is_num_user_missing (*vars, score)))
 603             {
 604               all_sysmis = 1;
 605               break;
 606             }
 607         }
 608     }
 609
 610   for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
 611     {
 612       double input = case_num (c, z->src_idx);
 613       double *output = &case_data_rw (c, z->dst_idx)->f;
 614
 615       if (z->mean == SYSMIS || z->std_dev == SYSMIS
 616           || all_sysmis || input == SYSMIS
 617           || (!t->include_user_missing
 618               && var_is_num_user_missing (z->v, input)))
 619         *output = SYSMIS;
 620       else
 621         *output = (input - z->mean) / z->std_dev;
 622     }
 623   return TRNS_CONTINUE;
 624 }
 625
 626 /* Frees a descriptives_trns struct. */
 627 static bool
 628 descriptives_trns_free (void *trns_)
 629 {
 630   struct dsc_trns *t = trns_;
 631
 632   free (t->z_scores);
 633   assert((t->missing_type != DSC_LISTWISE) ^ (t->vars != NULL));
 634   free (t->vars);
 635   return true;
 636 }
 637
 638 /* Sets up a transformation to calculate Z scores. */
 639 static void
 640 setup_z_trns (struct dsc_proc *dsc, struct dataset *ds)
 641 {
 642   struct dsc_trns *t;
 643   size_t cnt, i;
 644
 645   for (cnt = i = 0; i < dsc->var_cnt; i++)
 646     if (dsc->vars[i].z_name[0] != '\0')
 647       cnt++;
 648
 649   t = xmalloc (sizeof *t);
 650   t->z_scores = xnmalloc (cnt, sizeof *t->z_scores);
 651   t->z_score_cnt = cnt;
 652   t->missing_type = dsc->missing_type;
 653   t->include_user_missing = dsc->include_user_missing;
 654   if ( t->missing_type == DSC_LISTWISE )
 655     {
 656       t->var_cnt = dsc->var_cnt;
 657       t->vars = xnmalloc (t->var_cnt, sizeof *t->vars);
 658       for (i = 0; i < t->var_cnt; i++)
 659         t->vars[i] = dsc->vars[i].v;
 660     }
 661   else
 662     {
 663       t->var_cnt = 0;
 664       t->vars = NULL;
 665     }
 666
 667   for (cnt = i = 0; i < dsc->var_cnt; i++)
 668     {
 669       struct dsc_var *dv = &dsc->vars[i];
 670       if (dv->z_name[0] != '\0')
 671         {
 672           struct dsc_z_score *z;
 673           struct variable *dst_var;
 674
 675           dst_var = dict_create_var_assert (dataset_dict (ds), dv->z_name, 0);
 676           var_set_label (dst_var, xasprintf (_("Z-score of %s"),
 677                                              var_to_string (dv->v)));
 678
 679           z = &t->z_scores[cnt++];
 680           z->src_idx = dv->v->fv;
 681           z->dst_idx = dst_var->fv;
 682           z->mean = dv->stats[DSC_MEAN];
 683           z->std_dev = dv->stats[DSC_STDDEV];
 684           z->v = dv->v;
 685         }
 686     }
 687
 688   add_transformation (ds,
 689                       descriptives_trns_proc, descriptives_trns_free, t);
 690 }
 691 \f
 692 /* Statistical calculation. */
 693
 694 static bool listwise_missing (struct dsc_proc *dsc, const struct ccase *c);
 695
 696 /* Calculates and displays descriptive statistics for the cases
 697    in CF. */
 698 static bool
 699 calc_descriptives (const struct ccase *first,
 700                    const struct casefile *cf, void *dsc_,
 701                    const struct dataset *ds)
 702 {
 703   struct dsc_proc *dsc = dsc_;
 704   struct casereader *reader;
 705   struct ccase c;
 706   size_t i;
 707
 708   output_split_file_values (ds, first);
 709
 710   for (i = 0; i < dsc->var_cnt; i++)
 711     {
 712       struct dsc_var *dv = &dsc->vars[i];
 713
 714       dv->valid = dv->missing = 0.0;
 715       if (dv->moments != NULL)
 716         moments_clear (dv->moments);
 717       dv->min = DBL_MAX;
 718       dv->max = -DBL_MAX;
 719     }
 720   dsc->missing_listwise = 0.;
 721   dsc->valid = 0.;
 722
 723   /* First pass to handle most of the work. */
 724   for (reader = casefile_get_reader (cf, NULL);
 725        casereader_read (reader, &c);
 726        case_destroy (&c))
 727     {
 728       double weight = dict_get_case_weight (dataset_dict (ds), &c, &dsc->bad_warn);
 729       if (weight <= 0.0)
 730         continue;
 731
 732       /* Check for missing values. */
 733       if (listwise_missing (dsc, &c))
 734         {
 735           dsc->missing_listwise += weight;
 736           if (dsc->missing_type == DSC_LISTWISE)
 737             continue;
 738         }
 739       dsc->valid += weight;
 740
 741       for (i = 0; i < dsc->var_cnt; i++)
 742         {
 743           struct dsc_var *dv = &dsc->vars[i];
 744           double x = case_num (&c, dv->v->fv);
 745
 746           if (dsc->missing_type != DSC_LISTWISE
 747               && (x == SYSMIS
 748                   || (!dsc->include_user_missing
 749                       && var_is_num_user_missing (dv->v, x))))
 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   casereader_destroy (reader);
 765
 766   /* Second pass for higher-order moments. */
 767   if (dsc->max_moment > MOMENT_MEAN)
 768     {
 769       for (reader = casefile_get_reader (cf, NULL);
 770            casereader_read (reader, &c);
 771            case_destroy (&c))
 772         {
 773           double weight = dict_get_case_weight (dataset_dict (ds), &c,
 774                                                 &dsc->bad_warn);
 775           if (weight <= 0.0)
 776             continue;
 777
 778           /* Check for missing values. */
 779           if (listwise_missing (dsc, &c)
 780               && dsc->missing_type == DSC_LISTWISE)
 781             continue;
 782
 783           for (i = 0; i < dsc->var_cnt; i++)
 784             {
 785               struct dsc_var *dv = &dsc->vars[i];
 786               double x = case_num (&c, dv->v->fv);
 787
 788               if (dsc->missing_type != DSC_LISTWISE
 789                   && (x == SYSMIS
 790                       || (!dsc->include_user_missing
 791                           && var_is_num_user_missing (dv->v, x))))
 792                 continue;
 793
 794               if (dv->moments != NULL)
 795                 moments_pass_two (dv->moments, x, weight);
 796             }
 797         }
 798       casereader_destroy (reader);
 799     }
 800
 801   /* Calculate results. */
 802   for (i = 0; i < dsc->var_cnt; i++)
 803     {
 804       struct dsc_var *dv = &dsc->vars[i];
 805       double W;
 806       int j;
 807
 808       for (j = 0; j < DSC_N_STATS; j++)
 809         dv->stats[j] = SYSMIS;
 810
 811       dv->valid = W = dsc->valid - dv->missing;
 812
 813       if (dv->moments != NULL)
 814         moments_calculate (dv->moments, NULL,
 815                            &dv->stats[DSC_MEAN], &dv->stats[DSC_VARIANCE],
 816                            &dv->stats[DSC_SKEWNESS], &dv->stats[DSC_KURTOSIS]);
 817       if (dsc->calc_stats & (1ul << DSC_SEMEAN)
 818           && dv->stats[DSC_VARIANCE] != SYSMIS && W > 0.)
 819         dv->stats[DSC_SEMEAN] = sqrt (dv->stats[DSC_VARIANCE]) / sqrt (W);
 820       if (dsc->calc_stats & (1ul << DSC_STDDEV)
 821           && dv->stats[DSC_VARIANCE] != SYSMIS)
 822         dv->stats[DSC_STDDEV] = sqrt (dv->stats[DSC_VARIANCE]);
 823       if (dsc->calc_stats & (1ul << DSC_SEKURT))
 824         if (dv->stats[DSC_KURTOSIS] != SYSMIS)
 825             dv->stats[DSC_SEKURT] = calc_sekurt (W);
 826       if (dsc->calc_stats & (1ul << DSC_SESKEW)
 827           && dv->stats[DSC_SKEWNESS] != SYSMIS)
 828         dv->stats[DSC_SESKEW] = calc_seskew (W);
 829       dv->stats[DSC_RANGE] = ((dv->min == DBL_MAX || dv->max == -DBL_MAX)
 830                               ? SYSMIS : dv->max - dv->min);
 831       dv->stats[DSC_MIN] = dv->min == DBL_MAX ? SYSMIS : dv->min;
 832       dv->stats[DSC_MAX] = dv->max == -DBL_MAX ? SYSMIS : dv->max;
 833       if (dsc->calc_stats & (1ul << DSC_SUM))
 834         dv->stats[DSC_SUM] = W * dv->stats[DSC_MEAN];
 835     }
 836
 837   /* Output results. */
 838   display (dsc);
 839
 840   return true;
 841 }
 842
 843 /* Returns true if any of the descriptives variables in DSC's
 844    variable list have missing values in case C, false otherwise. */
 845 static bool
 846 listwise_missing (struct dsc_proc *dsc, const struct ccase *c)
 847 {
 848   size_t i;
 849
 850   for (i = 0; i < dsc->var_cnt; i++)
 851     {
 852       struct dsc_var *dv = &dsc->vars[i];
 853       double x = case_num (c, dv->v->fv);
 854
 855       if (x == SYSMIS
 856           || (!dsc->include_user_missing
 857               && var_is_num_user_missing (dv->v, x)))
 858         return true;
 859     }
 860   return false;
 861 }
 862 \f
 863 /* Statistical display. */
 864
 865 static algo_compare_func descriptives_compare_dsc_vars;
 866
 867 /* Displays a table of descriptive statistics for DSC. */
 868 static void
 869 display (struct dsc_proc *dsc)
 870 {
 871   size_t i;
 872   int nc;
 873   struct tab_table *t;
 874
 875   nc = 1 + (dsc->format == DSC_SERIAL ? 2 : 1);
 876   for (i = 0; i < DSC_N_STATS; i++)
 877     if (dsc->show_stats & (1ul << i))
 878       nc++;
 879
 880   if (dsc->sort_by_stat != DSC_NONE)
 881     sort (dsc->vars, dsc->var_cnt, sizeof *dsc->vars,
 882           descriptives_compare_dsc_vars, dsc);
 883
 884   t = tab_create (nc, dsc->var_cnt + 1, 0);
 885   tab_headers (t, 1, 0, 1, 0);
 886   tab_box (t, TAL_1, TAL_1, -1, -1, 0, 0, nc - 1, dsc->var_cnt);
 887   tab_box (t, -1, -1, -1, TAL_1, 1, 0, nc - 1, dsc->var_cnt);
 888   tab_hline (t, TAL_2, 0, nc - 1, 1);
 889   tab_vline (t, TAL_2, 1, 0, dsc->var_cnt);
 890   tab_dim (t, tab_natural_dimensions);
 891
 892   nc = 0;
 893   tab_text (t, nc++, 0, TAB_LEFT | TAT_TITLE, _("Variable"));
 894   if (dsc->format == DSC_SERIAL)
 895     {
 896       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Valid N"));
 897       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Missing N"));
 898     }
 899   else
 900     tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, "N");
 901
 902   for (i = 0; i < DSC_N_STATS; i++)
 903     if (dsc->show_stats & (1ul << i))
 904       {
 905         const char *title = gettext (dsc_info[i].name);
 906         tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, title);
 907       }
 908
 909   for (i = 0; i < dsc->var_cnt; i++)
 910     {
 911       struct dsc_var *dv = &dsc->vars[i];
 912       size_t j;
 913
 914       nc = 0;
 915       tab_text (t, nc++, i + 1, TAB_LEFT, var_get_name (dv->v));
 916       tab_text (t, nc++, i + 1, TAT_PRINTF, "%g", dv->valid);
 917       if (dsc->format == DSC_SERIAL)
 918         tab_text (t, nc++, i + 1, TAT_PRINTF, "%g", dv->missing);
 919       for (j = 0; j < DSC_N_STATS; j++)
 920         if (dsc->show_stats & (1ul << j))
 921           tab_float (t, nc++, i + 1, TAB_NONE, dv->stats[j], 10, 3);
 922     }
 923
 924   tab_title (t, _("Valid cases = %g; cases with missing value(s) = %g."),
 925              dsc->valid, dsc->missing_listwise);
 926
 927   tab_submit (t);
 928 }
 929
 930 /* Compares `struct dsc_var's A and B according to the ordering
 931    specified by CMD. */
 932 static int
 933 descriptives_compare_dsc_vars (const void *a_, const void *b_, const void *dsc_)
 934 {
 935   const struct dsc_var *a = a_;
 936   const struct dsc_var *b = b_;
 937   const struct dsc_proc *dsc = dsc_;
 938
 939   int result;
 940
 941   if (dsc->sort_by_stat == DSC_NAME)
 942     result = strcasecmp (var_get_name (a->v), var_get_name (b->v));
 943   else
 944     {
 945       double as = a->stats[dsc->sort_by_stat];
 946       double bs = b->stats[dsc->sort_by_stat];
 947
 948       result = as < bs ? -1 : as > bs;
 949     }
 950
 951   if (!dsc->sort_ascending)
 952     result = -result;
 953
 954   return result;
 955 }