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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 1997-2000, 2009-2014 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 <float.h>
  20 #include <limits.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23
  24 #include "data/casegrouper.h"
  25 #include "data/casereader.h"
  26 #include "data/casewriter.h"
  27 #include "data/dataset.h"
  28 #include "data/dictionary.h"
  29 #include "data/transformations.h"
  30 #include "data/variable.h"
  31 #include "language/command.h"
  32 #include "language/dictionary/split-file.h"
  33 #include "language/lexer/lexer.h"
  34 #include "language/lexer/variable-parser.h"
  35 #include "libpspp/array.h"
  36 #include "libpspp/assertion.h"
  37 #include "libpspp/compiler.h"
  38 #include "libpspp/i18n.h"
  39 #include "libpspp/message.h"
  40 #include "math/moments.h"
  41 #include "output/pivot-table.h"
  42
  43 #include "gl/xalloc.h"
  44
  45 #include "gettext.h"
  46 #define _(msgid) gettext (msgid)
  47 #define N_(msgid) msgid
  48
  49 /* DESCRIPTIVES private data. */
  50
  51 struct dsc_proc;
  52
  53 /* Handling of missing values. */
  54 enum dsc_missing_type
  55   {
  56     DSC_VARIABLE,       /* Handle missing values on a per-variable basis. */
  57     DSC_LISTWISE        /* Discard entire case if any variable is missing. */
  58   };
  59
  60 /* Describes properties of a distribution for the purpose of
  61    calculating a Z-score. */
  62 struct dsc_z_score
  63   {
  64     const struct variable *src_var;   /* Variable on which z-score is based. */
  65     struct variable *z_var;     /* New z-score variable. */
  66     double mean;                /* Distribution mean. */
  67     double std_dev;             /* Distribution standard deviation. */
  68   };
  69
  70 /* DESCRIPTIVES transformation (for calculating Z-scores). */
  71 struct dsc_trns
  72   {
  73     struct dsc_z_score *z_scores; /* Array of Z-scores. */
  74     int n_z_scores;               /* Number of Z-scores. */
  75     const struct variable **vars;     /* Variables for listwise missing checks. */
  76     size_t n_vars;              /* Number of variables. */
  77     enum dsc_missing_type missing_type; /* Treatment of missing values. */
  78     enum mv_class exclude;      /* Classes of missing values to exclude. */
  79     const struct variable *filter;    /* Dictionary FILTER BY variable. */
  80     struct casereader *z_reader; /* Reader for count, mean, stddev. */
  81     casenumber count;            /* Number left in this SPLIT FILE group.*/
  82     bool ok;
  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     const struct variable *v;         /* Variable to calculate on. */
 132     char *z_name;                     /* 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 /* A DESCRIPTIVES procedure. */
 140 struct dsc_proc
 141   {
 142     /* Per-variable info. */
 143     struct dictionary *dict;    /* Dictionary. */
 144     struct dsc_var *vars;       /* Variables. */
 145     size_t n_vars;              /* 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
 151     /* Accumulated results. */
 152     double missing_listwise;    /* Sum of weights of cases missing listwise. */
 153     double valid;               /* Sum of weights of valid cases. */
 154     bool bad_warn;               /* Warn if bad weight found. */
 155     enum dsc_statistic sort_by_stat; /* Statistic to sort by; -1: name. */
 156     int sort_ascending;         /* !0: ascending order; 0: descending. */
 157     unsigned long show_stats;   /* Statistics to display. */
 158     unsigned long calc_stats;   /* Statistics to calculate. */
 159     enum moment max_moment;     /* Highest moment needed for stats. */
 160
 161     /* Z scores. */
 162     struct casewriter *z_writer; /* Mean and stddev per SPLIT FILE group. */
 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 *n_zs);
 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 n_vars = 0;
 193   int save_z_scores = 0;
 194   int n_zs = 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->dict = dict;
 204   dsc->vars = NULL;
 205   dsc->n_vars = 0;
 206   dsc->missing_type = DSC_VARIABLE;
 207   dsc->exclude = MV_ANY;
 208   dsc->missing_listwise = 0.;
 209   dsc->valid = 0.;
 210   dsc->bad_warn = 1;
 211   dsc->sort_by_stat = DSC_NONE;
 212   dsc->sort_ascending = 1;
 213   dsc->show_stats = dsc->calc_stats = DEFAULT_STATS;
 214   dsc->z_writer = NULL;
 215
 216   /* Parse DESCRIPTIVES. */
 217   int z_ofs = 0;
 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         {
 241           save_z_scores = 1;
 242           z_ofs = lex_ofs (lexer) - 1;
 243         }
 244       else if (lex_match_id (lexer, "FORMAT"))
 245         {
 246           lex_match (lexer, T_EQUALS);
 247           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 248             {
 249               if (lex_match_id (lexer, "LABELS")
 250                   || lex_match_id (lexer, "NOLABELS")
 251                   || lex_match_id (lexer, "INDEX")
 252                   || lex_match_id (lexer, "NOINDEX")
 253                   || lex_match_id (lexer, "LINE")
 254                   || lex_match_id (lexer, "SERIAL"))
 255                 {
 256                   /* Ignore. */
 257                 }
 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                 {
 278                   enum dsc_statistic s = match_statistic (lexer);
 279                   if (s == DSC_NONE)
 280                     {
 281                       lex_error (lexer, NULL);
 282                       goto error;
 283                     }
 284                   dsc->show_stats |= 1ul << s;
 285                 }
 286               lex_match (lexer, T_COMMA);
 287             }
 288           if (dsc->show_stats == 0)
 289             dsc->show_stats = DEFAULT_STATS;
 290         }
 291       else if (lex_match_id (lexer, "SORT"))
 292         {
 293           lex_match (lexer, T_EQUALS);
 294           if (lex_match_id (lexer, "NAME"))
 295             dsc->sort_by_stat = DSC_NAME;
 296           else
 297             {
 298               dsc->sort_by_stat = match_statistic (lexer);
 299               if (dsc->sort_by_stat == DSC_NONE)
 300                 dsc->sort_by_stat = DSC_MEAN;
 301             }
 302           if (lex_match (lexer, T_LPAREN))
 303             {
 304               if (lex_match_id (lexer, "A"))
 305                 dsc->sort_ascending = 1;
 306               else if (lex_match_id (lexer, "D"))
 307                 dsc->sort_ascending = 0;
 308               else
 309                 lex_error (lexer, NULL);
 310               if (! lex_force_match (lexer, T_RPAREN))
 311                 goto error;
 312             }
 313         }
 314       else if (n_vars == 0)
 315         {
 316           if (lex_next_token (lexer, 1) == T_EQUALS)
 317             {
 318               lex_match_id (lexer, "VARIABLES");
 319               lex_match (lexer, T_EQUALS);
 320             }
 321
 322           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 323             {
 324               int i;
 325
 326               if (!parse_variables_const (lexer, dict, &vars, &n_vars,
 327                                     PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
 328                 goto error;
 329
 330               dsc->vars = xnrealloc ((void *)dsc->vars, n_vars, sizeof *dsc->vars);
 331               for (i = dsc->n_vars; i < n_vars; i++)
 332                 {
 333                   struct dsc_var *dv = &dsc->vars[i];
 334                   dv->v = vars[i];
 335                   dv->z_name = NULL;
 336                   dv->moments = NULL;
 337                 }
 338               dsc->n_vars = n_vars;
 339
 340               if (lex_match (lexer, T_LPAREN))
 341                 {
 342                   if (!lex_force_id (lexer))
 343                     goto error;
 344                   z_ofs = lex_ofs (lexer);
 345                   if (try_name (dict, dsc, lex_tokcstr (lexer)))
 346                     {
 347                       struct dsc_var *dsc_var = &dsc->vars[dsc->n_vars - 1];
 348                       dsc_var->z_name = xstrdup (lex_tokcstr (lexer));
 349                       n_zs++;
 350                     }
 351                   else
 352                     lex_error (lexer, _("Z-score variable name %s would be "
 353                                         "a duplicate variable name."),
 354                                lex_tokcstr (lexer));
 355                   lex_get (lexer);
 356                   if (!lex_force_match (lexer, T_RPAREN))
 357                     goto error;
 358                 }
 359             }
 360         }
 361       else
 362         {
 363           lex_error (lexer, NULL);
 364           goto error;
 365         }
 366
 367       lex_match (lexer, T_SLASH);
 368     }
 369   if (n_vars == 0)
 370     {
 371       msg (SE, _("No variables specified."));
 372       goto error;
 373     }
 374
 375   /* Construct z-score varnames, show translation table. */
 376   if (n_zs || save_z_scores)
 377     {
 378       struct caseproto *proto;
 379
 380       if (save_z_scores)
 381         {
 382           int n_gens = 0;
 383
 384           for (i = 0; i < dsc->n_vars; i++)
 385             {
 386               struct dsc_var *dsc_var = &dsc->vars[i];
 387               if (dsc_var->z_name == NULL)
 388                 {
 389                   const char *name = var_get_name (dsc_var->v);
 390                   dsc_var->z_name = generate_z_varname (dict, dsc, name,
 391                                                         &n_gens);
 392                   if (dsc_var->z_name == NULL)
 393                     goto error;
 394
 395                   n_zs++;
 396                 }
 397             }
 398         }
 399
 400       /* It would be better to handle Z scores correctly (however we define
 401          that) when TEMPORARY is in effect, but in the meantime this at least
 402          prevents a use-after-free error.  See bug #38786.  */
 403       if (proc_make_temporary_transformations_permanent (ds))
 404         lex_ofs_msg (lexer, SW, z_ofs, z_ofs,
 405                      _("DESCRIPTIVES with Z scores ignores TEMPORARY.  "
 406                        "Temporary transformations will be made permanent."));
 407
 408       proto = caseproto_create ();
 409       for (i = 0; i < 1 + 2 * n_zs; i++)
 410         proto = caseproto_add_width (proto, 0);
 411       dsc->z_writer = autopaging_writer_create (proto);
 412       caseproto_unref (proto);
 413
 414       dump_z_table (dsc);
 415     }
 416
 417   /* Figure out statistics to display. */
 418   if (dsc->show_stats & (1ul << DSC_SKEWNESS))
 419     dsc->show_stats |= 1ul << DSC_SESKEW;
 420   if (dsc->show_stats & (1ul << DSC_KURTOSIS))
 421     dsc->show_stats |= 1ul << DSC_SEKURT;
 422
 423   /* Figure out which statistics to calculate. */
 424   dsc->calc_stats = dsc->show_stats;
 425   if (n_zs > 0)
 426     dsc->calc_stats |= (1ul << DSC_MEAN) | (1ul << DSC_STDDEV);
 427   if (dsc->sort_by_stat >= 0)
 428     dsc->calc_stats |= 1ul << dsc->sort_by_stat;
 429   if (dsc->show_stats & (1ul << DSC_SESKEW))
 430     dsc->calc_stats |= 1ul << DSC_SKEWNESS;
 431   if (dsc->show_stats & (1ul << DSC_SEKURT))
 432     dsc->calc_stats |= 1ul << DSC_KURTOSIS;
 433
 434   /* Figure out maximum moment needed and allocate moments for
 435      the variables. */
 436   dsc->max_moment = MOMENT_NONE;
 437   for (i = 0; i < DSC_N_STATS; i++)
 438     if (dsc->calc_stats & (1ul << i) && dsc_info[i].moment > dsc->max_moment)
 439       dsc->max_moment = dsc_info[i].moment;
 440   if (dsc->max_moment != MOMENT_NONE)
 441     for (i = 0; i < dsc->n_vars; i++)
 442       dsc->vars[i].moments = moments_create (dsc->max_moment);
 443
 444   /* Data pass. */
 445   grouper = casegrouper_create_splits (proc_open_filtering (ds, false), dict);
 446   while (casegrouper_get_next_group (grouper, &group))
 447     calc_descriptives (dsc, group, ds);
 448   ok = casegrouper_destroy (grouper);
 449   ok = proc_commit (ds) && ok;
 450
 451   /* Z-scoring! */
 452   if (ok && n_zs)
 453     setup_z_trns (dsc, ds);
 454
 455   /* Done. */
 456   free (vars);
 457   free_dsc_proc (dsc);
 458   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
 459
 460  error:
 461   free (vars);
 462   free_dsc_proc (dsc);
 463   return CMD_FAILURE;
 464 }
 465
 466 /* Returns the statistic named by the current token and skips past the token.
 467    Returns DSC_NONE if no statistic is given (e.g., subcommand with no
 468    specifiers). Emits an error if the current token ID does not name a
 469    statistic. */
 470 static enum dsc_statistic
 471 match_statistic (struct lexer *lexer)
 472 {
 473   if (lex_token (lexer) == T_ID)
 474     {
 475       for (enum dsc_statistic stat = 0; stat < DSC_N_STATS; stat++)
 476         if (lex_match_id (lexer, dsc_info[stat].identifier))
 477           return stat;
 478
 479       const char *stat_names[DSC_N_STATS];
 480       for (enum dsc_statistic stat = 0; stat < DSC_N_STATS; stat++)
 481         stat_names[stat] = dsc_info[stat].identifier;
 482       lex_error_expecting_array (lexer, stat_names,
 483                                  sizeof stat_names / sizeof *stat_names);
 484       lex_get (lexer);
 485     }
 486
 487   return DSC_NONE;
 488 }
 489
 490 /* Frees DSC. */
 491 static void
 492 free_dsc_proc (struct dsc_proc *dsc)
 493 {
 494   size_t i;
 495
 496   if (dsc == NULL)
 497     return;
 498
 499   for (i = 0; i < dsc->n_vars; i++)
 500     {
 501       struct dsc_var *dsc_var = &dsc->vars[i];
 502       free (dsc_var->z_name);
 503       moments_destroy (dsc_var->moments);
 504     }
 505   casewriter_destroy (dsc->z_writer);
 506   free (dsc->vars);
 507   free (dsc);
 508 }
 509 \f
 510 /* Z scores. */
 511
 512 /* Returns false if NAME is a duplicate of any existing variable name or
 513    of any previously-declared z-var name; otherwise returns true. */
 514 static bool
 515 try_name (const struct dictionary *dict, struct dsc_proc *dsc,
 516           const char *name)
 517 {
 518   size_t i;
 519
 520   if (dict_lookup_var (dict, name) != NULL)
 521     return false;
 522   for (i = 0; i < dsc->n_vars; i++)
 523     {
 524       struct dsc_var *dsc_var = &dsc->vars[i];
 525       if (dsc_var->z_name != NULL && !utf8_strcasecmp (dsc_var->z_name, name))
 526         return false;
 527     }
 528   return true;
 529 }
 530
 531 /* Generates a name for a Z-score variable based on a variable
 532    named VAR_NAME, given that *Z_CNT generated variable names are
 533    known to already exist.  If successful, returns the new name
 534    as a dynamically allocated string.  On failure, returns NULL. */
 535 static char *
 536 generate_z_varname (const struct dictionary *dict, struct dsc_proc *dsc,
 537                     const char *var_name, int *n_zs)
 538 {
 539   char *z_name, *trunc_name;
 540
 541   /* Try a name based on the original variable name. */
 542   z_name = xasprintf ("Z%s", var_name);
 543   trunc_name = utf8_encoding_trunc (z_name, dict_get_encoding (dict),
 544                                     ID_MAX_LEN);
 545   free (z_name);
 546   if (try_name (dict, dsc, trunc_name))
 547     return trunc_name;
 548   free (trunc_name);
 549
 550   /* Generate a synthetic name. */
 551   for (;;)
 552     {
 553       char name[16];
 554
 555       (*n_zs)++;
 556
 557       if (*n_zs <= 99)
 558         sprintf (name, "ZSC%03d", *n_zs);
 559       else if (*n_zs <= 108)
 560         sprintf (name, "STDZ%02d", *n_zs - 99);
 561       else if (*n_zs <= 117)
 562         sprintf (name, "ZZZZ%02d", *n_zs - 108);
 563       else if (*n_zs <= 126)
 564         sprintf (name, "ZQZQ%02d", *n_zs - 117);
 565       else
 566         {
 567           msg (SE, _("Ran out of generic names for Z-score variables.  "
 568                      "There are only 126 generic names: ZSC001-ZSC0999, "
 569                      "STDZ01-STDZ09, ZZZZ01-ZZZZ09, ZQZQ01-ZQZQ09."));
 570           return NULL;
 571         }
 572
 573       if (try_name (dict, dsc, name))
 574         return xstrdup (name);
 575     }
 576   NOT_REACHED();
 577 }
 578
 579 /* Outputs a table describing the mapping between source
 580    variables and Z-score variables. */
 581 static void
 582 dump_z_table (struct dsc_proc *dsc)
 583 {
 584   struct pivot_table *table = pivot_table_create (
 585     N_("Mapping of Variables to Z-scores"));
 586
 587   pivot_dimension_create (table, PIVOT_AXIS_COLUMN, N_("Names"),
 588                           N_("Source"), N_("Target"));
 589
 590   struct pivot_dimension *names = pivot_dimension_create (
 591     table, PIVOT_AXIS_ROW, N_("Variables"));
 592   names->hide_all_labels = true;
 593
 594   for (size_t i = 0; i < dsc->n_vars; i++)
 595     if (dsc->vars[i].z_name != NULL)
 596       {
 597         int row = pivot_category_create_leaf (names->root,
 598                                               pivot_value_new_number (i));
 599
 600         pivot_table_put2 (table, 0, row,
 601                           pivot_value_new_variable (dsc->vars[i].v));
 602         pivot_table_put2 (table, 1, row,
 603                           pivot_value_new_user_text (dsc->vars[i].z_name, -1));
 604       }
 605
 606   pivot_table_submit (table);
 607 }
 608
 609 static void
 610 descriptives_set_all_sysmis_zscores (const struct dsc_trns *t, struct ccase *c)
 611 {
 612   const struct dsc_z_score *z;
 613
 614   for (z = t->z_scores; z < t->z_scores + t->n_z_scores; z++)
 615     *case_num_rw (c, z->z_var) = SYSMIS;
 616 }
 617
 618 /* Transformation function to calculate Z-scores. Will return SYSMIS if any of
 619    the following are true: 1) mean or standard deviation is SYSMIS 2) score is
 620    SYSMIS 3) score is user missing and they were not included in the original
 621    analyis. 4) any of the variables in the original analysis were missing
 622    (either system or user-missing values that weren't included).
 623 */
 624 static enum trns_result
 625 descriptives_trns_proc (void *trns_, struct ccase **c,
 626                         casenumber case_idx UNUSED)
 627 {
 628   struct dsc_trns *t = trns_;
 629   struct dsc_z_score *z;
 630   const struct variable **vars;
 631
 632   *c = case_unshare (*c);
 633
 634   if (t->filter)
 635     {
 636       double f = case_num (*c, t->filter);
 637       if (f == 0.0 || var_is_num_missing (t->filter, f))
 638         {
 639           descriptives_set_all_sysmis_zscores (t, *c);
 640           return TRNS_CONTINUE;
 641         }
 642     }
 643
 644   if (t->count <= 0)
 645     {
 646       struct ccase *z_case;
 647
 648       z_case = casereader_read (t->z_reader);
 649       if (z_case)
 650         {
 651           size_t z_idx = 0;
 652
 653           t->count = case_num_idx (z_case, z_idx++);
 654           for (z = t->z_scores; z < t->z_scores + t->n_z_scores; z++)
 655             {
 656               z->mean = case_num_idx (z_case, z_idx++);
 657               z->std_dev = case_num_idx (z_case, z_idx++);
 658             }
 659           case_unref (z_case);
 660         }
 661       else
 662         {
 663           if (t->ok)
 664             {
 665               msg (SE,  _("Internal error processing Z scores.  "
 666                           "Please report this to %s."),
 667                    PACKAGE_BUGREPORT);
 668               t->ok = false;
 669             }
 670           descriptives_set_all_sysmis_zscores (t, *c);
 671           return TRNS_CONTINUE;
 672         }
 673     }
 674   t->count--;
 675
 676   if (t->missing_type == DSC_LISTWISE)
 677     {
 678       assert(t->vars);
 679       for (vars = t->vars; vars < t->vars + t->n_vars; vars++)
 680         {
 681           double score = case_num (*c, *vars);
 682           if (var_is_num_missing (*vars, score) & t->exclude)
 683             {
 684               descriptives_set_all_sysmis_zscores (t, *c);
 685               return TRNS_CONTINUE;
 686             }
 687         }
 688     }
 689
 690   for (z = t->z_scores; z < t->z_scores + t->n_z_scores; z++)
 691     {
 692       double input = case_num (*c, z->src_var);
 693       double *output = case_num_rw (*c, z->z_var);
 694
 695       if (z->mean == SYSMIS || z->std_dev == SYSMIS
 696           || var_is_num_missing (z->src_var, input) & t->exclude)
 697         *output = SYSMIS;
 698       else
 699         *output = (input - z->mean) / z->std_dev;
 700     }
 701   return TRNS_CONTINUE;
 702 }
 703
 704 /* Frees a descriptives_trns struct. */
 705 static bool
 706 descriptives_trns_free (void *trns_)
 707 {
 708   struct dsc_trns *t = trns_;
 709   bool ok = t->ok && !casereader_error (t->z_reader);
 710
 711   free (t->z_scores);
 712   casereader_destroy (t->z_reader);
 713   assert((t->missing_type != DSC_LISTWISE) ^ (t->vars != NULL));
 714   free (t->vars);
 715   free (t);
 716
 717   return ok;
 718 }
 719
 720 static const struct trns_class descriptives_trns_class = {
 721   .name = "DESCRIPTIVES (Z scores)",
 722   .execute = descriptives_trns_proc,
 723   .destroy = descriptives_trns_free,
 724 };
 725
 726 /* Sets up a transformation to calculate Z scores. */
 727 static void
 728 setup_z_trns (struct dsc_proc *dsc, struct dataset *ds)
 729 {
 730   struct dsc_trns *t;
 731   size_t n, i;
 732
 733   for (n = i = 0; i < dsc->n_vars; i++)
 734     if (dsc->vars[i].z_name != NULL)
 735       n++;
 736
 737   t = xmalloc (sizeof *t);
 738   t->z_scores = xnmalloc (n, sizeof *t->z_scores);
 739   t->n_z_scores = n;
 740   t->missing_type = dsc->missing_type;
 741   t->exclude = dsc->exclude;
 742   if (t->missing_type == DSC_LISTWISE)
 743     {
 744       t->n_vars = dsc->n_vars;
 745       t->vars = xnmalloc (t->n_vars, sizeof *t->vars);
 746       for (i = 0; i < t->n_vars; i++)
 747         t->vars[i] = dsc->vars[i].v;
 748     }
 749   else
 750     {
 751       t->n_vars = 0;
 752       t->vars = NULL;
 753     }
 754   t->filter = dict_get_filter (dataset_dict (ds));
 755   t->z_reader = casewriter_make_reader (dsc->z_writer);
 756   t->count = 0;
 757   t->ok = true;
 758   dsc->z_writer = NULL;
 759
 760   for (n = i = 0; i < dsc->n_vars; i++)
 761     {
 762       struct dsc_var *dv = &dsc->vars[i];
 763       if (dv->z_name != NULL)
 764         {
 765           struct dsc_z_score *z;
 766           struct variable *dst_var;
 767           char *label;
 768
 769           dst_var = dict_create_var_assert (dataset_dict (ds), dv->z_name, 0);
 770
 771           label = xasprintf (_("Z-score of %s"),var_to_string (dv->v));
 772           var_set_label (dst_var, label);
 773           free (label);
 774
 775           z = &t->z_scores[n++];
 776           z->src_var = dv->v;
 777           z->z_var = dst_var;
 778         }
 779     }
 780
 781   add_transformation (ds, &descriptives_trns_class, t);
 782 }
 783 \f
 784 /* Statistical calculation. */
 785
 786 static bool listwise_missing (struct dsc_proc *dsc, const struct ccase *c);
 787
 788 /* Calculates and displays descriptive statistics for the cases
 789    in CF. */
 790 static void
 791 calc_descriptives (struct dsc_proc *dsc, struct casereader *group,
 792                    struct dataset *ds)
 793 {
 794   const struct variable *filter = dict_get_filter (dataset_dict (ds));
 795   struct casereader *pass1, *pass2;
 796   casenumber count;
 797   struct ccase *c;
 798   size_t z_idx;
 799   size_t i;
 800
 801   c = casereader_peek (group, 0);
 802   if (c == NULL)
 803     {
 804       casereader_destroy (group);
 805       return;
 806     }
 807   output_split_file_values (ds, c);
 808   case_unref (c);
 809
 810   group = casereader_create_filter_weight (group, dataset_dict (ds),
 811                                            NULL, NULL);
 812
 813   pass1 = group;
 814   pass2 = dsc->max_moment <= MOMENT_MEAN ? NULL : casereader_clone (pass1);
 815
 816   for (i = 0; i < dsc->n_vars; i++)
 817     {
 818       struct dsc_var *dv = &dsc->vars[i];
 819
 820       dv->valid = dv->missing = 0.0;
 821       if (dv->moments != NULL)
 822         moments_clear (dv->moments);
 823       dv->min = DBL_MAX;
 824       dv->max = -DBL_MAX;
 825     }
 826   dsc->missing_listwise = 0.;
 827   dsc->valid = 0.;
 828
 829   /* First pass to handle most of the work. */
 830   count = 0;
 831   for (; (c = casereader_read (pass1)) != NULL; case_unref (c))
 832     {
 833       double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
 834
 835       if (filter)
 836         {
 837           double f = case_num (c, filter);
 838           if (f == 0.0 || var_is_num_missing (filter, f))
 839             continue;
 840         }
 841
 842       /* Check for missing values. */
 843       if (listwise_missing (dsc, c))
 844         {
 845           dsc->missing_listwise += weight;
 846           if (dsc->missing_type == DSC_LISTWISE)
 847             continue;
 848         }
 849       dsc->valid += weight;
 850
 851       for (i = 0; i < dsc->n_vars; i++)
 852         {
 853           struct dsc_var *dv = &dsc->vars[i];
 854           double x = case_num (c, dv->v);
 855
 856           if (var_is_num_missing (dv->v, x) & dsc->exclude)
 857             {
 858               dv->missing += weight;
 859               continue;
 860             }
 861
 862           if (dv->moments != NULL)
 863             moments_pass_one (dv->moments, x, weight);
 864
 865           if (x < dv->min)
 866             dv->min = x;
 867           if (x > dv->max)
 868             dv->max = x;
 869         }
 870
 871       count++;
 872     }
 873   if (!casereader_destroy (pass1))
 874     {
 875       casereader_destroy (pass2);
 876       return;
 877     }
 878
 879   /* Second pass for higher-order moments. */
 880   if (dsc->max_moment > MOMENT_MEAN)
 881     {
 882       for (; (c = casereader_read (pass2)) != NULL; case_unref (c))
 883         {
 884           double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
 885
 886           if (filter)
 887             {
 888               double f = case_num (c, filter);
 889               if (f == 0.0 || var_is_num_missing (filter, f))
 890                 continue;
 891             }
 892
 893           /* Check for missing values. */
 894           if (dsc->missing_type == DSC_LISTWISE && listwise_missing (dsc, c))
 895             continue;
 896
 897           for (i = 0; i < dsc->n_vars; i++)
 898             {
 899               struct dsc_var *dv = &dsc->vars[i];
 900               double x = case_num (c, dv->v);
 901
 902               if (var_is_num_missing (dv->v, x) & dsc->exclude)
 903                 continue;
 904
 905               if (dv->moments != NULL)
 906                 moments_pass_two (dv->moments, x, weight);
 907             }
 908         }
 909       if (!casereader_destroy (pass2))
 910         return;
 911     }
 912
 913   /* Calculate results. */
 914   if (dsc->z_writer && count > 0)
 915     {
 916       c = case_create (casewriter_get_proto (dsc->z_writer));
 917       z_idx = 0;
 918       *case_num_rw_idx (c, z_idx++) = count;
 919     }
 920   else
 921     c = NULL;
 922
 923   for (i = 0; i < dsc->n_vars; i++)
 924     {
 925       struct dsc_var *dv = &dsc->vars[i];
 926       double W;
 927       int j;
 928
 929       for (j = 0; j < DSC_N_STATS; j++)
 930         dv->stats[j] = SYSMIS;
 931
 932       dv->valid = W = dsc->valid - dv->missing;
 933
 934       if (dv->moments != NULL)
 935         moments_calculate (dv->moments, NULL,
 936                            &dv->stats[DSC_MEAN], &dv->stats[DSC_VARIANCE],
 937                            &dv->stats[DSC_SKEWNESS], &dv->stats[DSC_KURTOSIS]);
 938       if (dsc->calc_stats & (1ul << DSC_SEMEAN)
 939           && dv->stats[DSC_VARIANCE] != SYSMIS && W > 0.)
 940         dv->stats[DSC_SEMEAN] = sqrt (dv->stats[DSC_VARIANCE]) / sqrt (W);
 941       if (dsc->calc_stats & (1ul << DSC_STDDEV)
 942           && dv->stats[DSC_VARIANCE] != SYSMIS)
 943         dv->stats[DSC_STDDEV] = sqrt (dv->stats[DSC_VARIANCE]);
 944       if (dsc->calc_stats & (1ul << DSC_SEKURT))
 945         if (dv->stats[DSC_KURTOSIS] != SYSMIS)
 946             dv->stats[DSC_SEKURT] = calc_sekurt (W);
 947       if (dsc->calc_stats & (1ul << DSC_SESKEW)
 948           && dv->stats[DSC_SKEWNESS] != SYSMIS)
 949         dv->stats[DSC_SESKEW] = calc_seskew (W);
 950       dv->stats[DSC_RANGE] = ((dv->min == DBL_MAX || dv->max == -DBL_MAX)
 951                               ? SYSMIS : dv->max - dv->min);
 952       dv->stats[DSC_MIN] = dv->min == DBL_MAX ? SYSMIS : dv->min;
 953       dv->stats[DSC_MAX] = dv->max == -DBL_MAX ? SYSMIS : dv->max;
 954       if (dsc->calc_stats & (1ul << DSC_SUM))
 955         dv->stats[DSC_SUM] = W * dv->stats[DSC_MEAN];
 956
 957       if (dv->z_name && c != NULL)
 958         {
 959           *case_num_rw_idx (c, z_idx++) = dv->stats[DSC_MEAN];
 960           *case_num_rw_idx (c, z_idx++) = dv->stats[DSC_STDDEV];
 961         }
 962     }
 963
 964   if (c != NULL)
 965     casewriter_write (dsc->z_writer, c);
 966
 967   /* Output results. */
 968   display (dsc);
 969 }
 970
 971 /* Returns true if any of the descriptives variables in DSC's
 972    variable list have missing values in case C, false otherwise. */
 973 static bool
 974 listwise_missing (struct dsc_proc *dsc, const struct ccase *c)
 975 {
 976   size_t i;
 977
 978   for (i = 0; i < dsc->n_vars; i++)
 979     {
 980       struct dsc_var *dv = &dsc->vars[i];
 981       double x = case_num (c, dv->v);
 982
 983       if (var_is_num_missing (dv->v, x) & dsc->exclude)
 984         return true;
 985     }
 986   return false;
 987 }
 988 \f
 989 /* Statistical display. */
 990
 991 static algo_compare_func descriptives_compare_dsc_vars;
 992
 993 /* Displays a table of descriptive statistics for DSC. */
 994 static void
 995 display (struct dsc_proc *dsc)
 996 {
 997   struct pivot_table *table = pivot_table_create (
 998     N_("Descriptive Statistics"));
 999   pivot_table_set_weight_var (table, dict_get_weight (dsc->dict));
1000
1001   struct pivot_dimension *statistics = pivot_dimension_create (
1002     table, PIVOT_AXIS_COLUMN, N_("Statistics"));
1003   pivot_category_create_leaf_rc (
1004     statistics->root, pivot_value_new_text (N_("N")), PIVOT_RC_COUNT);
1005   for (int i = 0; i < DSC_N_STATS; i++)
1006     if (dsc->show_stats & (1ul << i))
1007       pivot_category_create_leaf (statistics->root,
1008                                   pivot_value_new_text (dsc_info[i].name));
1009
1010   if (dsc->sort_by_stat != DSC_NONE)
1011     sort (dsc->vars, dsc->n_vars, sizeof *dsc->vars,
1012           descriptives_compare_dsc_vars, dsc);
1013
1014   struct pivot_dimension *variables = pivot_dimension_create (
1015     table, PIVOT_AXIS_ROW, N_("Variable"));
1016   for (size_t i = 0; i < dsc->n_vars; i++)
1017     {
1018       const struct dsc_var *dv = &dsc->vars[i];
1019
1020       int row = pivot_category_create_leaf (variables->root,
1021                                             pivot_value_new_variable (dv->v));
1022
1023       int column = 0;
1024       pivot_table_put2 (table, column++, row,
1025                         pivot_value_new_number (dv->valid));
1026
1027       for (int j = 0; j < DSC_N_STATS; j++)
1028         if (dsc->show_stats & (1ul << j))
1029           {
1030             union value v = { .f = dv->stats[j] };
1031             struct pivot_value *pv = (j == DSC_MIN || j == DSC_MAX
1032                                       ? pivot_value_new_var_value (dv->v, &v)
1033                                       : pivot_value_new_number (dv->stats[j]));
1034             pivot_table_put2 (table, column++, row, pv);
1035           }
1036     }
1037
1038   int row = pivot_category_create_leaves (
1039     variables->root, N_("Valid N (listwise)"), N_("Missing N (listwise)"));
1040   pivot_table_put2 (table, 0, row, pivot_value_new_number (dsc->valid));
1041   pivot_table_put2 (table, 0, row + 1,
1042                     pivot_value_new_number (dsc->missing_listwise));
1043   pivot_table_submit (table);
1044 }
1045
1046 /* Compares `struct dsc_var's A and B according to the ordering
1047    specified by CMD. */
1048 static int
1049 descriptives_compare_dsc_vars (const void *a_, const void *b_, const void *dsc_)
1050 {
1051   const struct dsc_var *a = a_;
1052   const struct dsc_var *b = b_;
1053   const struct dsc_proc *dsc = dsc_;
1054
1055   int result;
1056
1057   if (dsc->sort_by_stat == DSC_NAME)
1058     result = utf8_strcasecmp (var_get_name (a->v), var_get_name (b->v));
1059   else
1060     {
1061       double as = a->stats[dsc->sort_by_stat];
1062       double bs = b->stats[dsc->sort_by_stat];
1063
1064       result = as < bs ? -1 : as > bs;
1065     }
1066
1067   if (!dsc->sort_ascending)
1068     result = -result;
1069
1070   return result;
1071 }