pintos-os.org Git - pspp/blob - src/language/stats/descriptives.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 1997-9, 2000, 2009, 2010, 2011 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/dataset.h"
  26 #include "data/dictionary.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/assertion.h"
  35 #include "libpspp/compiler.h"
  36 #include "libpspp/i18n.h"
  37 #include "libpspp/message.h"
  38 #include "math/moments.h"
  39 #include "output/tab.h"
  40
  41 #include "gl/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;                     /* 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 char *generate_z_varname (const struct dictionary *dict,
 173                                  struct dsc_proc *dsc,
 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) != T_ENDCMD)
 219     {
 220       if (lex_match_id (lexer, "MISSING"))
 221         {
 222           lex_match (lexer, T_EQUALS);
 223           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 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, T_COMMA);
 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, T_EQUALS);
 244           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 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, T_COMMA);
 264             }
 265         }
 266       else if (lex_match_id (lexer, "STATISTICS"))
 267         {
 268           lex_match (lexer, T_EQUALS);
 269           dsc->show_stats = 0;
 270           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 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, T_COMMA);
 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, T_EQUALS);
 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, T_LPAREN))
 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, T_RPAREN);
 303             }
 304         }
 305       else if (var_cnt == 0)
 306         {
 307           if (lex_next_token (lexer, 1) == T_EQUALS)
 308             {
 309               lex_match_id (lexer, "VARIABLES");
 310               lex_match (lexer, T_EQUALS);
 311             }
 312
 313           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 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 = NULL;
 327                   dv->moments = NULL;
 328                 }
 329               dsc->var_cnt = var_cnt;
 330
 331               if (lex_match (lexer, T_LPAREN))
 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_tokcstr (lexer)))
 339                     {
 340                       struct dsc_var *dsc_var = &dsc->vars[dsc->var_cnt - 1];
 341                       dsc_var->z_name = xstrdup (lex_tokcstr (lexer));
 342                       z_cnt++;
 343                     }
 344                   else
 345                     msg (SE, _("Z-score variable name %s would be"
 346                                " a duplicate variable name."), lex_tokcstr (lexer));
 347                   lex_get (lexer);
 348                   if (!lex_force_match (lexer, T_RPAREN))
 349                     goto error;
 350                 }
 351             }
 352         }
 353       else
 354         {
 355           lex_error (lexer, NULL);
 356           goto error;
 357         }
 358
 359       lex_match (lexer, T_SLASH);
 360     }
 361   if (var_cnt == 0)
 362     {
 363       msg (SE, _("No variables specified."));
 364       goto error;
 365     }
 366
 367   /* Construct z-score varnames, show translation table. */
 368   if (z_cnt || save_z_scores)
 369     {
 370       if (save_z_scores)
 371         {
 372           int gen_cnt = 0;
 373
 374           for (i = 0; i < dsc->var_cnt; i++)
 375             {
 376               struct dsc_var *dsc_var = &dsc->vars[i];
 377               if (dsc_var->z_name == NULL)
 378                 {
 379                   const char *name = var_get_name (dsc_var->v);
 380                   dsc_var->z_name = generate_z_varname (dict, dsc, name,
 381                                                         &gen_cnt);
 382                   if (dsc_var->z_name == NULL)
 383                     goto error;
 384
 385                   z_cnt++;
 386                 }
 387             }
 388         }
 389       dump_z_table (dsc);
 390     }
 391
 392   /* Figure out statistics to display. */
 393   if (dsc->show_stats & (1ul << DSC_SKEWNESS))
 394     dsc->show_stats |= 1ul << DSC_SESKEW;
 395   if (dsc->show_stats & (1ul << DSC_KURTOSIS))
 396     dsc->show_stats |= 1ul << DSC_SEKURT;
 397
 398   /* Figure out which statistics to calculate. */
 399   dsc->calc_stats = dsc->show_stats;
 400   if (z_cnt > 0)
 401     dsc->calc_stats |= (1ul << DSC_MEAN) | (1ul << DSC_STDDEV);
 402   if (dsc->sort_by_stat >= 0)
 403     dsc->calc_stats |= 1ul << dsc->sort_by_stat;
 404   if (dsc->show_stats & (1ul << DSC_SESKEW))
 405     dsc->calc_stats |= 1ul << DSC_SKEWNESS;
 406   if (dsc->show_stats & (1ul << DSC_SEKURT))
 407     dsc->calc_stats |= 1ul << DSC_KURTOSIS;
 408
 409   /* Figure out maximum moment needed and allocate moments for
 410      the variables. */
 411   dsc->max_moment = MOMENT_NONE;
 412   for (i = 0; i < DSC_N_STATS; i++)
 413     if (dsc->calc_stats & (1ul << i) && dsc_info[i].moment > dsc->max_moment)
 414       dsc->max_moment = dsc_info[i].moment;
 415   if (dsc->max_moment != MOMENT_NONE)
 416     for (i = 0; i < dsc->var_cnt; i++)
 417       dsc->vars[i].moments = moments_create (dsc->max_moment);
 418
 419   /* Data pass. */
 420   grouper = casegrouper_create_splits (proc_open (ds), dict);
 421   while (casegrouper_get_next_group (grouper, &group))
 422     calc_descriptives (dsc, group, ds);
 423   ok = casegrouper_destroy (grouper);
 424   ok = proc_commit (ds) && ok;
 425
 426   /* Z-scoring! */
 427   if (ok && z_cnt)
 428     setup_z_trns (dsc, ds);
 429
 430   /* Done. */
 431   free (vars);
 432   free_dsc_proc (dsc);
 433   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
 434
 435  error:
 436   free (vars);
 437   free_dsc_proc (dsc);
 438   return CMD_FAILURE;
 439 }
 440
 441 /* Returns the statistic named by the current token and skips past the token.
 442    Returns DSC_NONE if no statistic is given (e.g., subcommand with no
 443    specifiers). Emits an error if the current token ID does not name a
 444    statistic. */
 445 static enum dsc_statistic
 446 match_statistic (struct lexer *lexer)
 447 {
 448   if (lex_token (lexer) == T_ID)
 449     {
 450       enum dsc_statistic stat;
 451
 452       for (stat = 0; stat < DSC_N_STATS; stat++)
 453         if (lex_match_id (lexer, dsc_info[stat].identifier))
 454           return stat;
 455
 456       lex_get (lexer);
 457       lex_error (lexer, _("expecting statistic name: reverting to default"));
 458     }
 459
 460   return DSC_NONE;
 461 }
 462
 463 /* Frees DSC. */
 464 static void
 465 free_dsc_proc (struct dsc_proc *dsc)
 466 {
 467   size_t i;
 468
 469   if (dsc == NULL)
 470     return;
 471
 472   for (i = 0; i < dsc->var_cnt; i++)
 473     {
 474       struct dsc_var *dsc_var = &dsc->vars[i];
 475       free (dsc_var->z_name);
 476       moments_destroy (dsc_var->moments);
 477     }
 478   free (dsc->vars);
 479   free (dsc);
 480 }
 481 \f
 482 /* Z scores. */
 483
 484 /* Returns false if NAME is a duplicate of any existing variable name or
 485    of any previously-declared z-var name; otherwise returns true. */
 486 static bool
 487 try_name (const struct dictionary *dict, struct dsc_proc *dsc,
 488           const char *name)
 489 {
 490   size_t i;
 491
 492   if (dict_lookup_var (dict, name) != NULL)
 493     return false;
 494   for (i = 0; i < dsc->var_cnt; i++)
 495     {
 496       struct dsc_var *dsc_var = &dsc->vars[i];
 497       if (dsc_var->z_name != NULL && !strcasecmp (dsc_var->z_name, name))
 498         return false;
 499     }
 500   return true;
 501 }
 502
 503 /* Generates a name for a Z-score variable based on a variable
 504    named VAR_NAME, given that *Z_CNT generated variable names are
 505    known to already exist.  If successful, returns the new name
 506    as a dynamically allocated string.  On failure, returns NULL. */
 507 static char *
 508 generate_z_varname (const struct dictionary *dict, struct dsc_proc *dsc,
 509                     const char *var_name, int *z_cnt)
 510 {
 511   char *z_name, *trunc_name;
 512
 513   /* Try a name based on the original variable name. */
 514   z_name = xasprintf ("Z%s", var_name);
 515   trunc_name = utf8_encoding_trunc (z_name, dict_get_encoding (dict),
 516                                     ID_MAX_LEN);
 517   free (z_name);
 518   if (try_name (dict, dsc, trunc_name))
 519     return trunc_name;
 520   free (trunc_name);
 521
 522   /* Generate a synthetic name. */
 523   for (;;)
 524     {
 525       char name[8];
 526
 527       (*z_cnt)++;
 528
 529       if (*z_cnt <= 99)
 530         sprintf (name, "ZSC%03d", *z_cnt);
 531       else if (*z_cnt <= 108)
 532         sprintf (name, "STDZ%02d", *z_cnt - 99);
 533       else if (*z_cnt <= 117)
 534         sprintf (name, "ZZZZ%02d", *z_cnt - 108);
 535       else if (*z_cnt <= 126)
 536         sprintf (name, "ZQZQ%02d", *z_cnt - 117);
 537       else
 538         {
 539           msg (SE, _("Ran out of generic names for Z-score variables.  "
 540                      "There are only 126 generic names: ZSC001-ZSC0999, "
 541                      "STDZ01-STDZ09, ZZZZ01-ZZZZ09, ZQZQ01-ZQZQ09."));
 542           return NULL;
 543         }
 544
 545       if (try_name (dict, dsc, name))
 546         return xstrdup (name);
 547     }
 548   NOT_REACHED();
 549 }
 550
 551 /* Outputs a table describing the mapping between source
 552    variables and Z-score variables. */
 553 static void
 554 dump_z_table (struct dsc_proc *dsc)
 555 {
 556   size_t cnt = 0;
 557   struct tab_table *t;
 558
 559   {
 560     size_t i;
 561
 562     for (i = 0; i < dsc->var_cnt; i++)
 563       if (dsc->vars[i].z_name != NULL)
 564         cnt++;
 565   }
 566
 567   t = tab_create (2, cnt + 1);
 568   tab_title (t, _("Mapping of variables to corresponding Z-scores."));
 569   tab_headers (t, 0, 0, 1, 0);
 570   tab_box (t, TAL_1, TAL_1, TAL_0, TAL_1, 0, 0, 1, cnt);
 571   tab_hline (t, TAL_2, 0, 1, 1);
 572   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
 573   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Target"));
 574
 575   {
 576     size_t i, y;
 577
 578     for (i = 0, y = 1; i < dsc->var_cnt; i++)
 579       if (dsc->vars[i].z_name != NULL)
 580         {
 581           tab_text (t, 0, y, TAB_LEFT, var_to_string (dsc->vars[i].v));
 582           tab_text (t, 1, y++, TAB_LEFT, dsc->vars[i].z_name);
 583         }
 584   }
 585
 586   tab_submit (t);
 587 }
 588
 589 /* Transformation function to calculate Z-scores. Will return SYSMIS if any of
 590    the following are true: 1) mean or standard deviation is SYSMIS 2) score is
 591    SYSMIS 3) score is user missing and they were not included in the original
 592    analyis. 4) any of the variables in the original analysis were missing
 593    (either system or user-missing values that weren't included).
 594 */
 595 static int
 596 descriptives_trns_proc (void *trns_, struct ccase **c,
 597                         casenumber case_idx UNUSED)
 598 {
 599   struct dsc_trns *t = trns_;
 600   struct dsc_z_score *z;
 601   const struct variable **vars;
 602   int all_sysmis = 0;
 603
 604   if (t->missing_type == DSC_LISTWISE)
 605     {
 606       assert(t->vars);
 607       for (vars = t->vars; vars < t->vars + t->var_cnt; vars++)
 608         {
 609           double score = case_num (*c, *vars);
 610           if (var_is_num_missing (*vars, score, t->exclude))
 611             {
 612               all_sysmis = 1;
 613               break;
 614             }
 615         }
 616     }
 617
 618   *c = case_unshare (*c);
 619   for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
 620     {
 621       double input = case_num (*c, z->src_var);
 622       double *output = &case_data_rw (*c, z->z_var)->f;
 623
 624       if (z->mean == SYSMIS || z->std_dev == SYSMIS || all_sysmis
 625           || var_is_num_missing (z->src_var, input, t->exclude))
 626         *output = SYSMIS;
 627       else
 628         *output = (input - z->mean) / z->std_dev;
 629     }
 630   return TRNS_CONTINUE;
 631 }
 632
 633 /* Frees a descriptives_trns struct. */
 634 static bool
 635 descriptives_trns_free (void *trns_)
 636 {
 637   struct dsc_trns *t = trns_;
 638
 639   free (t->z_scores);
 640   assert((t->missing_type != DSC_LISTWISE) ^ (t->vars != NULL));
 641   free (t->vars);
 642   return true;
 643 }
 644
 645 /* Sets up a transformation to calculate Z scores. */
 646 static void
 647 setup_z_trns (struct dsc_proc *dsc, struct dataset *ds)
 648 {
 649   struct dsc_trns *t;
 650   size_t cnt, i;
 651
 652   for (cnt = i = 0; i < dsc->var_cnt; i++)
 653     if (dsc->vars[i].z_name != NULL)
 654       cnt++;
 655
 656   t = xmalloc (sizeof *t);
 657   t->z_scores = xnmalloc (cnt, sizeof *t->z_scores);
 658   t->z_score_cnt = cnt;
 659   t->missing_type = dsc->missing_type;
 660   t->exclude = dsc->exclude;
 661   if ( t->missing_type == DSC_LISTWISE )
 662     {
 663       t->var_cnt = dsc->var_cnt;
 664       t->vars = xnmalloc (t->var_cnt, sizeof *t->vars);
 665       for (i = 0; i < t->var_cnt; i++)
 666         t->vars[i] = dsc->vars[i].v;
 667     }
 668   else
 669     {
 670       t->var_cnt = 0;
 671       t->vars = NULL;
 672     }
 673
 674   for (cnt = i = 0; i < dsc->var_cnt; i++)
 675     {
 676       struct dsc_var *dv = &dsc->vars[i];
 677       if (dv->z_name != NULL)
 678         {
 679           struct dsc_z_score *z;
 680           struct variable *dst_var;
 681
 682           dst_var = dict_create_var_assert (dataset_dict (ds), dv->z_name, 0);
 683           var_set_label (dst_var,
 684                          xasprintf (_("Z-score of %s"),var_to_string (dv->v)),
 685                          false);
 686
 687           z = &t->z_scores[cnt++];
 688           z->src_var = dv->v;
 689           z->z_var = dst_var;
 690           z->mean = dv->stats[DSC_MEAN];
 691           z->std_dev = dv->stats[DSC_STDDEV];
 692         }
 693     }
 694
 695   add_transformation (ds,
 696                       descriptives_trns_proc, descriptives_trns_free, t);
 697 }
 698 \f
 699 /* Statistical calculation. */
 700
 701 static bool listwise_missing (struct dsc_proc *dsc, const struct ccase *c);
 702
 703 /* Calculates and displays descriptive statistics for the cases
 704    in CF. */
 705 static void
 706 calc_descriptives (struct dsc_proc *dsc, struct casereader *group,
 707                    struct dataset *ds)
 708 {
 709   struct casereader *pass1, *pass2;
 710   struct ccase *c;
 711   size_t i;
 712
 713   c = casereader_peek (group, 0);
 714   if (c == NULL)
 715     {
 716       casereader_destroy (group);
 717       return;
 718     }
 719   output_split_file_values (ds, c);
 720   case_unref (c);
 721
 722   group = casereader_create_filter_weight (group, dataset_dict (ds),
 723                                            NULL, NULL);
 724
 725   pass1 = group;
 726   pass2 = dsc->max_moment <= MOMENT_MEAN ? NULL : casereader_clone (pass1);
 727
 728   for (i = 0; i < dsc->var_cnt; i++)
 729     {
 730       struct dsc_var *dv = &dsc->vars[i];
 731
 732       dv->valid = dv->missing = 0.0;
 733       if (dv->moments != NULL)
 734         moments_clear (dv->moments);
 735       dv->min = DBL_MAX;
 736       dv->max = -DBL_MAX;
 737     }
 738   dsc->missing_listwise = 0.;
 739   dsc->valid = 0.;
 740
 741   /* First pass to handle most of the work. */
 742   for (; (c = casereader_read (pass1)) != NULL; case_unref (c))
 743     {
 744       double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
 745
 746       /* Check for missing values. */
 747       if (listwise_missing (dsc, c))
 748         {
 749           dsc->missing_listwise += weight;
 750           if (dsc->missing_type == DSC_LISTWISE)
 751             continue;
 752         }
 753       dsc->valid += weight;
 754
 755       for (i = 0; i < dsc->var_cnt; i++)
 756         {
 757           struct dsc_var *dv = &dsc->vars[i];
 758           double x = case_num (c, dv->v);
 759
 760           if (var_is_num_missing (dv->v, x, dsc->exclude))
 761             {
 762               dv->missing += weight;
 763               continue;
 764             }
 765
 766           if (dv->moments != NULL)
 767             moments_pass_one (dv->moments, x, weight);
 768
 769           if (x < dv->min)
 770             dv->min = x;
 771           if (x > dv->max)
 772             dv->max = x;
 773         }
 774     }
 775   if (!casereader_destroy (pass1))
 776     {
 777       casereader_destroy (pass2);
 778       return;
 779     }
 780
 781   /* Second pass for higher-order moments. */
 782   if (dsc->max_moment > MOMENT_MEAN)
 783     {
 784       for (; (c = casereader_read (pass2)) != NULL; case_unref (c))
 785         {
 786           double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
 787
 788           /* Check for missing values. */
 789           if (dsc->missing_type == DSC_LISTWISE && listwise_missing (dsc, c))
 790             continue;
 791
 792           for (i = 0; i < dsc->var_cnt; i++)
 793             {
 794               struct dsc_var *dv = &dsc->vars[i];
 795               double x = case_num (c, dv->v);
 796
 797               if (var_is_num_missing (dv->v, x, dsc->exclude))
 798                 continue;
 799
 800               if (dv->moments != NULL)
 801                 moments_pass_two (dv->moments, x, weight);
 802             }
 803         }
 804       if (!casereader_destroy (pass2))
 805         return;
 806     }
 807
 808   /* Calculate results. */
 809   for (i = 0; i < dsc->var_cnt; i++)
 810     {
 811       struct dsc_var *dv = &dsc->vars[i];
 812       double W;
 813       int j;
 814
 815       for (j = 0; j < DSC_N_STATS; j++)
 816         dv->stats[j] = SYSMIS;
 817
 818       dv->valid = W = dsc->valid - dv->missing;
 819
 820       if (dv->moments != NULL)
 821         moments_calculate (dv->moments, NULL,
 822                            &dv->stats[DSC_MEAN], &dv->stats[DSC_VARIANCE],
 823                            &dv->stats[DSC_SKEWNESS], &dv->stats[DSC_KURTOSIS]);
 824       if (dsc->calc_stats & (1ul << DSC_SEMEAN)
 825           && dv->stats[DSC_VARIANCE] != SYSMIS && W > 0.)
 826         dv->stats[DSC_SEMEAN] = sqrt (dv->stats[DSC_VARIANCE]) / sqrt (W);
 827       if (dsc->calc_stats & (1ul << DSC_STDDEV)
 828           && dv->stats[DSC_VARIANCE] != SYSMIS)
 829         dv->stats[DSC_STDDEV] = sqrt (dv->stats[DSC_VARIANCE]);
 830       if (dsc->calc_stats & (1ul << DSC_SEKURT))
 831         if (dv->stats[DSC_KURTOSIS] != SYSMIS)
 832             dv->stats[DSC_SEKURT] = calc_sekurt (W);
 833       if (dsc->calc_stats & (1ul << DSC_SESKEW)
 834           && dv->stats[DSC_SKEWNESS] != SYSMIS)
 835         dv->stats[DSC_SESKEW] = calc_seskew (W);
 836       dv->stats[DSC_RANGE] = ((dv->min == DBL_MAX || dv->max == -DBL_MAX)
 837                               ? SYSMIS : dv->max - dv->min);
 838       dv->stats[DSC_MIN] = dv->min == DBL_MAX ? SYSMIS : dv->min;
 839       dv->stats[DSC_MAX] = dv->max == -DBL_MAX ? SYSMIS : dv->max;
 840       if (dsc->calc_stats & (1ul << DSC_SUM))
 841         dv->stats[DSC_SUM] = W * dv->stats[DSC_MEAN];
 842     }
 843
 844   /* Output results. */
 845   display (dsc);
 846 }
 847
 848 /* Returns true if any of the descriptives variables in DSC's
 849    variable list have missing values in case C, false otherwise. */
 850 static bool
 851 listwise_missing (struct dsc_proc *dsc, const struct ccase *c)
 852 {
 853   size_t i;
 854
 855   for (i = 0; i < dsc->var_cnt; i++)
 856     {
 857       struct dsc_var *dv = &dsc->vars[i];
 858       double x = case_num (c, dv->v);
 859
 860       if (var_is_num_missing (dv->v, x, dsc->exclude))
 861         return true;
 862     }
 863   return false;
 864 }
 865 \f
 866 /* Statistical display. */
 867
 868 static algo_compare_func descriptives_compare_dsc_vars;
 869
 870 /* Displays a table of descriptive statistics for DSC. */
 871 static void
 872 display (struct dsc_proc *dsc)
 873 {
 874   size_t i;
 875   int nc;
 876   struct tab_table *t;
 877
 878   nc = 1 + (dsc->format == DSC_SERIAL ? 2 : 1);
 879   for (i = 0; i < DSC_N_STATS; i++)
 880     if (dsc->show_stats & (1ul << i))
 881       nc++;
 882
 883   if (dsc->sort_by_stat != DSC_NONE)
 884     sort (dsc->vars, dsc->var_cnt, sizeof *dsc->vars,
 885           descriptives_compare_dsc_vars, dsc);
 886
 887   t = tab_create (nc, dsc->var_cnt + 1);
 888   tab_headers (t, 1, 0, 1, 0);
 889   tab_box (t, TAL_1, TAL_1, -1, -1, 0, 0, nc - 1, dsc->var_cnt);
 890   tab_box (t, -1, -1, -1, TAL_1, 1, 0, nc - 1, dsc->var_cnt);
 891   tab_hline (t, TAL_2, 0, nc - 1, 1);
 892   tab_vline (t, TAL_2, 1, 0, dsc->var_cnt);
 893
 894   nc = 0;
 895   tab_text (t, nc++, 0, TAB_LEFT | TAT_TITLE, _("Variable"));
 896   if (dsc->format == DSC_SERIAL)
 897     {
 898       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Valid N"));
 899       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Missing N"));
 900     }
 901   else
 902     tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, "N");
 903
 904   for (i = 0; i < DSC_N_STATS; i++)
 905     if (dsc->show_stats & (1ul << i))
 906       {
 907         const char *title = gettext (dsc_info[i].name);
 908         tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, title);
 909       }
 910
 911   for (i = 0; i < dsc->var_cnt; i++)
 912     {
 913       struct dsc_var *dv = &dsc->vars[i];
 914       size_t j;
 915
 916       nc = 0;
 917       tab_text (t, nc++, i + 1, TAB_LEFT, var_to_string (dv->v));
 918       tab_text_format (t, nc++, i + 1, 0, "%g", dv->valid);
 919       if (dsc->format == DSC_SERIAL)
 920         tab_text_format (t, nc++, i + 1, 0, "%g", dv->missing);
 921
 922       for (j = 0; j < DSC_N_STATS; j++)
 923         if (dsc->show_stats & (1ul << j))
 924           tab_double (t, nc++, i + 1, TAB_NONE, dv->stats[j], NULL);
 925     }
 926
 927   tab_title (t, _("Valid cases = %g; cases with missing value(s) = %g."),
 928              dsc->valid, dsc->missing_listwise);
 929
 930   tab_submit (t);
 931 }
 932
 933 /* Compares `struct dsc_var's A and B according to the ordering
 934    specified by CMD. */
 935 static int
 936 descriptives_compare_dsc_vars (const void *a_, const void *b_, const void *dsc_)
 937 {
 938   const struct dsc_var *a = a_;
 939   const struct dsc_var *b = b_;
 940   const struct dsc_proc *dsc = dsc_;
 941
 942   int result;
 943
 944   if (dsc->sort_by_stat == DSC_NAME)
 945     result = strcasecmp (var_get_name (a->v), var_get_name (b->v));
 946   else
 947     {
 948       double as = a->stats[dsc->sort_by_stat];
 949       double bs = b->stats[dsc->sort_by_stat];
 950
 951       result = as < bs ? -1 : as > bs;
 952     }
 953
 954   if (!dsc->sort_ascending)
 955     result = -result;
 956
 957   return result;
 958 }