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

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 1997-2000, 2009-2013 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/casewriter.h"
  26 #include "data/dataset.h"
  27 #include "data/dictionary.h"
  28 #include "data/transformations.h"
  29 #include "data/variable.h"
  30 #include "language/command.h"
  31 #include "language/dictionary/split-file.h"
  32 #include "language/lexer/lexer.h"
  33 #include "language/lexer/variable-parser.h"
  34 #include "libpspp/array.h"
  35 #include "libpspp/assertion.h"
  36 #include "libpspp/compiler.h"
  37 #include "libpspp/i18n.h"
  38 #include "libpspp/message.h"
  39 #include "math/moments.h"
  40 #include "output/tab.h"
  41
  42 #include "gl/xalloc.h"
  43
  44 #include "gettext.h"
  45 #define _(msgid) gettext (msgid)
  46 #define N_(msgid) msgid
  47
  48 /* DESCRIPTIVES private data. */
  49
  50 struct dsc_proc;
  51
  52 /* Handling of missing values. */
  53 enum dsc_missing_type
  54   {
  55     DSC_VARIABLE,       /* Handle missing values on a per-variable basis. */
  56     DSC_LISTWISE        /* Discard entire case if any variable is missing. */
  57   };
  58
  59 /* Describes properties of a distribution for the purpose of
  60    calculating a Z-score. */
  61 struct dsc_z_score
  62   {
  63     const struct variable *src_var;   /* Variable on which z-score is based. */
  64     struct variable *z_var;     /* New z-score variable. */
  65     double mean;                /* Distribution mean. */
  66     double std_dev;             /* Distribution standard deviation. */
  67   };
  68
  69 /* DESCRIPTIVES transformation (for calculating Z-scores). */
  70 struct dsc_trns
  71   {
  72     struct dsc_z_score *z_scores; /* Array of Z-scores. */
  73     int z_score_cnt;            /* Number of Z-scores. */
  74     const struct variable **vars;     /* Variables for listwise missing checks. */
  75     size_t var_cnt;             /* Number of variables. */
  76     enum dsc_missing_type missing_type; /* Treatment of missing values. */
  77     enum mv_class exclude;      /* Classes of missing values to exclude. */
  78     struct casereader *z_reader; /* Reader for count, mean, stddev. */
  79     casenumber count;            /* Number left in this SPLIT FILE group.*/
  80     bool ok;
  81   };
  82
  83 /* Statistics.  Used as bit indexes, so must be 32 or fewer. */
  84 enum dsc_statistic
  85   {
  86     DSC_MEAN = 0, DSC_SEMEAN, DSC_STDDEV, DSC_VARIANCE, DSC_KURTOSIS,
  87     DSC_SEKURT, DSC_SKEWNESS, DSC_SESKEW, DSC_RANGE, DSC_MIN,
  88     DSC_MAX, DSC_SUM, DSC_N_STATS,
  89
  90     /* Only valid as sort criteria. */
  91     DSC_NAME = -2,              /* Sort by name. */
  92     DSC_NONE = -1               /* Unsorted. */
  93   };
  94
  95 /* Describes one statistic. */
  96 struct dsc_statistic_info
  97   {
  98     const char *identifier;     /* Identifier. */
  99     const char *name;           /* Full name. */
 100     enum moment moment;         /* Highest moment needed to calculate. */
 101   };
 102
 103 /* Table of statistics, indexed by DSC_*. */
 104 static const struct dsc_statistic_info dsc_info[DSC_N_STATS] =
 105   {
 106     {"MEAN", N_("Mean"), MOMENT_MEAN},
 107     {"SEMEAN", N_("S.E. Mean"), MOMENT_VARIANCE},
 108     {"STDDEV", N_("Std Dev"), MOMENT_VARIANCE},
 109     {"VARIANCE", N_("Variance"), MOMENT_VARIANCE},
 110     {"KURTOSIS", N_("Kurtosis"), MOMENT_KURTOSIS},
 111     {"SEKURTOSIS", N_("S.E. Kurt"), MOMENT_NONE},
 112     {"SKEWNESS", N_("Skewness"), MOMENT_SKEWNESS},
 113     {"SESKEWNESS", N_("S.E. Skew"), MOMENT_NONE},
 114     {"RANGE", N_("Range"), MOMENT_NONE},
 115     {"MINIMUM", N_("Minimum"), MOMENT_NONE},
 116     {"MAXIMUM", N_("Maximum"), MOMENT_NONE},
 117     {"SUM", N_("Sum"), MOMENT_MEAN},
 118   };
 119
 120 /* Statistics calculated by default if none are explicitly
 121    requested. */
 122 #define DEFAULT_STATS                                                   \
 123         ((1ul << DSC_MEAN) | (1ul << DSC_STDDEV) | (1ul << DSC_MIN)     \
 124          | (1ul << DSC_MAX))
 125
 126 /* A variable specified on DESCRIPTIVES. */
 127 struct dsc_var
 128   {
 129     const struct variable *v;         /* Variable to calculate on. */
 130     char *z_name;                     /* Name for z-score variable. */
 131     double valid, missing;      /* Valid, missing counts. */
 132     struct moments *moments;    /* Moments. */
 133     double min, max;            /* Maximum and mimimum values. */
 134     double stats[DSC_N_STATS];  /* All the stats' values. */
 135   };
 136
 137 /* Output format. */
 138 enum dsc_format
 139   {
 140     DSC_LINE,           /* Abbreviated format. */
 141     DSC_SERIAL          /* Long format. */
 142   };
 143
 144 /* A DESCRIPTIVES procedure. */
 145 struct dsc_proc
 146   {
 147     /* Per-variable info. */
 148     struct dsc_var *vars;       /* Variables. */
 149     size_t var_cnt;             /* Number of variables. */
 150
 151     /* User options. */
 152     enum dsc_missing_type missing_type; /* Treatment of missing values. */
 153     enum mv_class exclude;      /* Classes of missing values to exclude. */
 154     int show_var_labels;        /* Nonzero to show variable labels. */
 155     int show_index;             /* Nonzero to show variable index. */
 156     enum dsc_format format;     /* Output format. */
 157
 158     /* Accumulated results. */
 159     double missing_listwise;    /* Sum of weights of cases missing listwise. */
 160     double valid;               /* Sum of weights of valid cases. */
 161     bool bad_warn;               /* Warn if bad weight found. */
 162     enum dsc_statistic sort_by_stat; /* Statistic to sort by; -1: name. */
 163     int sort_ascending;         /* !0: ascending order; 0: descending. */
 164     unsigned long show_stats;   /* Statistics to display. */
 165     unsigned long calc_stats;   /* Statistics to calculate. */
 166     enum moment max_moment;     /* Highest moment needed for stats. */
 167
 168     /* Z scores. */
 169     struct casewriter *z_writer; /* Mean and stddev per SPLIT FILE group. */
 170   };
 171
 172 /* Parsing. */
 173 static enum dsc_statistic match_statistic (struct lexer *);
 174 static void free_dsc_proc (struct dsc_proc *);
 175
 176 /* Z-score functions. */
 177 static bool try_name (const struct dictionary *dict,
 178                       struct dsc_proc *dsc, const char *name);
 179 static char *generate_z_varname (const struct dictionary *dict,
 180                                  struct dsc_proc *dsc,
 181                                  const char *name, int *z_cnt);
 182 static void dump_z_table (struct dsc_proc *);
 183 static void setup_z_trns (struct dsc_proc *, struct dataset *);
 184
 185 /* Procedure execution functions. */
 186 static void calc_descriptives (struct dsc_proc *, struct casereader *,
 187                                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   const struct variable **vars = NULL;
 199   size_t var_cnt = 0;
 200   int save_z_scores = 0;
 201   int z_cnt = 0;
 202   size_t i;
 203   bool ok;
 204
 205   struct casegrouper *grouper;
 206   struct casereader *group;
 207
 208   /* Create and initialize dsc. */
 209   dsc = xmalloc (sizeof *dsc);
 210   dsc->vars = NULL;
 211   dsc->var_cnt = 0;
 212   dsc->missing_type = DSC_VARIABLE;
 213   dsc->exclude = MV_ANY;
 214   dsc->show_var_labels = 1;
 215   dsc->show_index = 0;
 216   dsc->format = DSC_LINE;
 217   dsc->missing_listwise = 0.;
 218   dsc->valid = 0.;
 219   dsc->bad_warn = 1;
 220   dsc->sort_by_stat = DSC_NONE;
 221   dsc->sort_ascending = 1;
 222   dsc->show_stats = dsc->calc_stats = DEFAULT_STATS;
 223   dsc->z_writer = NULL;
 224
 225   /* Parse DESCRIPTIVES. */
 226   while (lex_token (lexer) != T_ENDCMD)
 227     {
 228       if (lex_match_id (lexer, "MISSING"))
 229         {
 230           lex_match (lexer, T_EQUALS);
 231           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 232             {
 233               if (lex_match_id (lexer, "VARIABLE"))
 234                 dsc->missing_type = DSC_VARIABLE;
 235               else if (lex_match_id (lexer, "LISTWISE"))
 236                 dsc->missing_type = DSC_LISTWISE;
 237               else if (lex_match_id (lexer, "INCLUDE"))
 238                 dsc->exclude = MV_SYSTEM;
 239               else
 240                 {
 241                   lex_error (lexer, NULL);
 242                   goto error;
 243                 }
 244               lex_match (lexer, T_COMMA);
 245             }
 246         }
 247       else if (lex_match_id (lexer, "SAVE"))
 248         save_z_scores = 1;
 249       else if (lex_match_id (lexer, "FORMAT"))
 250         {
 251           lex_match (lexer, T_EQUALS);
 252           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 253             {
 254               if (lex_match_id (lexer, "LABELS"))
 255                 dsc->show_var_labels = 1;
 256               else if (lex_match_id (lexer, "NOLABELS"))
 257                 dsc->show_var_labels = 0;
 258               else if (lex_match_id (lexer, "INDEX"))
 259                 dsc->show_index = 1;
 260               else if (lex_match_id (lexer, "NOINDEX"))
 261                 dsc->show_index = 0;
 262               else if (lex_match_id (lexer, "LINE"))
 263                 dsc->format = DSC_LINE;
 264               else if (lex_match_id (lexer, "SERIAL"))
 265                 dsc->format = DSC_SERIAL;
 266               else
 267                 {
 268                   lex_error (lexer, NULL);
 269                   goto error;
 270                 }
 271               lex_match (lexer, T_COMMA);
 272             }
 273         }
 274       else if (lex_match_id (lexer, "STATISTICS"))
 275         {
 276           lex_match (lexer, T_EQUALS);
 277           dsc->show_stats = 0;
 278           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 279             {
 280               if (lex_match (lexer, T_ALL))
 281                 dsc->show_stats |= (1ul << DSC_N_STATS) - 1;
 282               else if (lex_match_id (lexer, "DEFAULT"))
 283                 dsc->show_stats |= DEFAULT_STATS;
 284               else
 285                 dsc->show_stats |= 1ul << (match_statistic (lexer));
 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               lex_force_match (lexer, T_RPAREN);
 311             }
 312         }
 313       else if (var_cnt == 0)
 314         {
 315           if (lex_next_token (lexer, 1) == T_EQUALS)
 316             {
 317               lex_match_id (lexer, "VARIABLES");
 318               lex_match (lexer, T_EQUALS);
 319             }
 320
 321           while (lex_token (lexer) != T_ENDCMD && lex_token (lexer) != T_SLASH)
 322             {
 323               int i;
 324
 325               if (!parse_variables_const (lexer, dict, &vars, &var_cnt,
 326                                     PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
 327                 goto error;
 328
 329               dsc->vars = xnrealloc ((void *)dsc->vars, var_cnt, sizeof *dsc->vars);
 330               for (i = dsc->var_cnt; i < var_cnt; i++)
 331                 {
 332                   struct dsc_var *dv = &dsc->vars[i];
 333                   dv->v = vars[i];
 334                   dv->z_name = NULL;
 335                   dv->moments = NULL;
 336                 }
 337               dsc->var_cnt = var_cnt;
 338
 339               if (lex_match (lexer, T_LPAREN))
 340                 {
 341                   if (lex_token (lexer) != T_ID)
 342                     {
 343                       lex_error (lexer, NULL);
 344                       goto error;
 345                     }
 346                   if (try_name (dict, dsc, lex_tokcstr (lexer)))
 347                     {
 348                       struct dsc_var *dsc_var = &dsc->vars[dsc->var_cnt - 1];
 349                       dsc_var->z_name = xstrdup (lex_tokcstr (lexer));
 350                       z_cnt++;
 351                     }
 352                   else
 353                     msg (SE, _("Z-score variable name %s would be"
 354                                " a duplicate variable name."), 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 (var_cnt == 0)
 370     {
 371       msg (SE, _("No variables specified."));
 372       goto error;
 373     }
 374
 375   /* Construct z-score varnames, show translation table. */
 376   if (z_cnt || save_z_scores)
 377     {
 378       struct caseproto *proto;
 379
 380       if (save_z_scores)
 381         {
 382           int gen_cnt = 0;
 383
 384           for (i = 0; i < dsc->var_cnt; 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                                                         &gen_cnt);
 392                   if (dsc_var->z_name == NULL)
 393                     goto error;
 394
 395                   z_cnt++;
 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         msg (SW, _("DESCRIPTIVES with Z scores ignores TEMPORARY.  "
 405                    "Temporary transformations will be made permanent."));
 406
 407       proto = caseproto_create ();
 408       for (i = 0; i < 1 + 2 * z_cnt; i++)
 409         proto = caseproto_add_width (proto, 0);
 410       dsc->z_writer = autopaging_writer_create (proto);
 411       caseproto_unref (proto);
 412
 413       dump_z_table (dsc);
 414     }
 415
 416   /* Figure out statistics to display. */
 417   if (dsc->show_stats & (1ul << DSC_SKEWNESS))
 418     dsc->show_stats |= 1ul << DSC_SESKEW;
 419   if (dsc->show_stats & (1ul << DSC_KURTOSIS))
 420     dsc->show_stats |= 1ul << DSC_SEKURT;
 421
 422   /* Figure out which statistics to calculate. */
 423   dsc->calc_stats = dsc->show_stats;
 424   if (z_cnt > 0)
 425     dsc->calc_stats |= (1ul << DSC_MEAN) | (1ul << DSC_STDDEV);
 426   if (dsc->sort_by_stat >= 0)
 427     dsc->calc_stats |= 1ul << dsc->sort_by_stat;
 428   if (dsc->show_stats & (1ul << DSC_SESKEW))
 429     dsc->calc_stats |= 1ul << DSC_SKEWNESS;
 430   if (dsc->show_stats & (1ul << DSC_SEKURT))
 431     dsc->calc_stats |= 1ul << DSC_KURTOSIS;
 432
 433   /* Figure out maximum moment needed and allocate moments for
 434      the variables. */
 435   dsc->max_moment = MOMENT_NONE;
 436   for (i = 0; i < DSC_N_STATS; i++)
 437     if (dsc->calc_stats & (1ul << i) && dsc_info[i].moment > dsc->max_moment)
 438       dsc->max_moment = dsc_info[i].moment;
 439   if (dsc->max_moment != MOMENT_NONE)
 440     for (i = 0; i < dsc->var_cnt; i++)
 441       dsc->vars[i].moments = moments_create (dsc->max_moment);
 442
 443   /* Data pass. */
 444   grouper = casegrouper_create_splits (proc_open_filtering (ds, z_cnt == 0),
 445                                        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 && z_cnt)
 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       enum dsc_statistic stat;
 476
 477       for (stat = 0; stat < DSC_N_STATS; stat++)
 478         if (lex_match_id (lexer, dsc_info[stat].identifier))
 479           return stat;
 480
 481       lex_get (lexer);
 482       lex_error (lexer, _("expecting statistic name: reverting to default"));
 483     }
 484
 485   return DSC_NONE;
 486 }
 487
 488 /* Frees DSC. */
 489 static void
 490 free_dsc_proc (struct dsc_proc *dsc)
 491 {
 492   size_t i;
 493
 494   if (dsc == NULL)
 495     return;
 496
 497   for (i = 0; i < dsc->var_cnt; i++)
 498     {
 499       struct dsc_var *dsc_var = &dsc->vars[i];
 500       free (dsc_var->z_name);
 501       moments_destroy (dsc_var->moments);
 502     }
 503   casewriter_destroy (dsc->z_writer);
 504   free (dsc->vars);
 505   free (dsc);
 506 }
 507 \f
 508 /* Z scores. */
 509
 510 /* Returns false if NAME is a duplicate of any existing variable name or
 511    of any previously-declared z-var name; otherwise returns true. */
 512 static bool
 513 try_name (const struct dictionary *dict, struct dsc_proc *dsc,
 514           const char *name)
 515 {
 516   size_t i;
 517
 518   if (dict_lookup_var (dict, name) != NULL)
 519     return false;
 520   for (i = 0; i < dsc->var_cnt; i++)
 521     {
 522       struct dsc_var *dsc_var = &dsc->vars[i];
 523       if (dsc_var->z_name != NULL && !utf8_strcasecmp (dsc_var->z_name, name))
 524         return false;
 525     }
 526   return true;
 527 }
 528
 529 /* Generates a name for a Z-score variable based on a variable
 530    named VAR_NAME, given that *Z_CNT generated variable names are
 531    known to already exist.  If successful, returns the new name
 532    as a dynamically allocated string.  On failure, returns NULL. */
 533 static char *
 534 generate_z_varname (const struct dictionary *dict, struct dsc_proc *dsc,
 535                     const char *var_name, int *z_cnt)
 536 {
 537   char *z_name, *trunc_name;
 538
 539   /* Try a name based on the original variable name. */
 540   z_name = xasprintf ("Z%s", var_name);
 541   trunc_name = utf8_encoding_trunc (z_name, dict_get_encoding (dict),
 542                                     ID_MAX_LEN);
 543   free (z_name);
 544   if (try_name (dict, dsc, trunc_name))
 545     return trunc_name;
 546   free (trunc_name);
 547
 548   /* Generate a synthetic name. */
 549   for (;;)
 550     {
 551       char name[8];
 552
 553       (*z_cnt)++;
 554
 555       if (*z_cnt <= 99)
 556         sprintf (name, "ZSC%03d", *z_cnt);
 557       else if (*z_cnt <= 108)
 558         sprintf (name, "STDZ%02d", *z_cnt - 99);
 559       else if (*z_cnt <= 117)
 560         sprintf (name, "ZZZZ%02d", *z_cnt - 108);
 561       else if (*z_cnt <= 126)
 562         sprintf (name, "ZQZQ%02d", *z_cnt - 117);
 563       else
 564         {
 565           msg (SE, _("Ran out of generic names for Z-score variables.  "
 566                      "There are only 126 generic names: ZSC001-ZSC0999, "
 567                      "STDZ01-STDZ09, ZZZZ01-ZZZZ09, ZQZQ01-ZQZQ09."));
 568           return NULL;
 569         }
 570
 571       if (try_name (dict, dsc, name))
 572         return xstrdup (name);
 573     }
 574   NOT_REACHED();
 575 }
 576
 577 /* Outputs a table describing the mapping between source
 578    variables and Z-score variables. */
 579 static void
 580 dump_z_table (struct dsc_proc *dsc)
 581 {
 582   size_t cnt = 0;
 583   struct tab_table *t;
 584
 585   {
 586     size_t i;
 587
 588     for (i = 0; i < dsc->var_cnt; i++)
 589       if (dsc->vars[i].z_name != NULL)
 590         cnt++;
 591   }
 592
 593   t = tab_create (2, cnt + 1);
 594   tab_title (t, _("Mapping of variables to corresponding Z-scores."));
 595   tab_headers (t, 0, 0, 1, 0);
 596   tab_box (t, TAL_1, TAL_1, TAL_0, TAL_1, 0, 0, 1, cnt);
 597   tab_hline (t, TAL_2, 0, 1, 1);
 598   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
 599   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Target"));
 600
 601   {
 602     size_t i, y;
 603
 604     for (i = 0, y = 1; i < dsc->var_cnt; i++)
 605       if (dsc->vars[i].z_name != NULL)
 606         {
 607           tab_text (t, 0, y, TAB_LEFT, var_to_string (dsc->vars[i].v));
 608           tab_text (t, 1, y++, TAB_LEFT, dsc->vars[i].z_name);
 609         }
 610   }
 611
 612   tab_submit (t);
 613 }
 614
 615 /* Transformation function to calculate Z-scores. Will return SYSMIS if any of
 616    the following are true: 1) mean or standard deviation is SYSMIS 2) score is
 617    SYSMIS 3) score is user missing and they were not included in the original
 618    analyis. 4) any of the variables in the original analysis were missing
 619    (either system or user-missing values that weren't included).
 620 */
 621 static int
 622 descriptives_trns_proc (void *trns_, struct ccase **c,
 623                         casenumber case_idx UNUSED)
 624 {
 625   struct dsc_trns *t = trns_;
 626   struct dsc_z_score *z;
 627   const struct variable **vars;
 628   int all_sysmis = 0;
 629
 630   if (t->count <= 0)
 631     {
 632       struct ccase *z_case;
 633
 634       z_case = casereader_read (t->z_reader);
 635       if (z_case)
 636         {
 637           size_t z_idx = 0;
 638
 639           t->count = case_num_idx (z_case, z_idx++);
 640           for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
 641             {
 642               z->mean = case_num_idx (z_case, z_idx++);
 643               z->std_dev = case_num_idx (z_case, z_idx++);
 644             }
 645           case_unref (z_case);
 646         }
 647       else
 648         {
 649           if (t->ok)
 650             {
 651               msg (SE, _("Internal error processing Z scores"));
 652               t->ok = false;
 653             }
 654           for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
 655             z->mean = z->std_dev = SYSMIS;
 656         }
 657     }
 658   t->count--;
 659
 660   if (t->missing_type == DSC_LISTWISE)
 661     {
 662       assert(t->vars);
 663       for (vars = t->vars; vars < t->vars + t->var_cnt; vars++)
 664         {
 665           double score = case_num (*c, *vars);
 666           if (var_is_num_missing (*vars, score, t->exclude))
 667             {
 668               all_sysmis = 1;
 669               break;
 670             }
 671         }
 672     }
 673
 674   *c = case_unshare (*c);
 675   for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
 676     {
 677       double input = case_num (*c, z->src_var);
 678       double *output = &case_data_rw (*c, z->z_var)->f;
 679
 680       if (z->mean == SYSMIS || z->std_dev == SYSMIS || all_sysmis
 681           || var_is_num_missing (z->src_var, input, t->exclude))
 682         *output = SYSMIS;
 683       else
 684         *output = (input - z->mean) / z->std_dev;
 685     }
 686   return TRNS_CONTINUE;
 687 }
 688
 689 /* Frees a descriptives_trns struct. */
 690 static bool
 691 descriptives_trns_free (void *trns_)
 692 {
 693   struct dsc_trns *t = trns_;
 694   bool ok = t->ok && !casereader_error (t->z_reader);
 695
 696   free (t->z_scores);
 697   casereader_destroy (t->z_reader);
 698   assert((t->missing_type != DSC_LISTWISE) ^ (t->vars != NULL));
 699   free (t->vars);
 700   free (t);
 701
 702   return ok;
 703 }
 704
 705 /* Sets up a transformation to calculate Z scores. */
 706 static void
 707 setup_z_trns (struct dsc_proc *dsc, struct dataset *ds)
 708 {
 709   struct dsc_trns *t;
 710   size_t cnt, i;
 711
 712   for (cnt = i = 0; i < dsc->var_cnt; i++)
 713     if (dsc->vars[i].z_name != NULL)
 714       cnt++;
 715
 716   t = xmalloc (sizeof *t);
 717   t->z_scores = xnmalloc (cnt, sizeof *t->z_scores);
 718   t->z_score_cnt = cnt;
 719   t->missing_type = dsc->missing_type;
 720   t->exclude = dsc->exclude;
 721   if ( t->missing_type == DSC_LISTWISE )
 722     {
 723       t->var_cnt = dsc->var_cnt;
 724       t->vars = xnmalloc (t->var_cnt, sizeof *t->vars);
 725       for (i = 0; i < t->var_cnt; i++)
 726         t->vars[i] = dsc->vars[i].v;
 727     }
 728   else
 729     {
 730       t->var_cnt = 0;
 731       t->vars = NULL;
 732     }
 733   t->z_reader = casewriter_make_reader (dsc->z_writer);
 734   t->count = 0;
 735   t->ok = true;
 736   dsc->z_writer = NULL;
 737
 738   for (cnt = i = 0; i < dsc->var_cnt; i++)
 739     {
 740       struct dsc_var *dv = &dsc->vars[i];
 741       if (dv->z_name != NULL)
 742         {
 743           struct dsc_z_score *z;
 744           struct variable *dst_var;
 745           char *label;
 746
 747           dst_var = dict_create_var_assert (dataset_dict (ds), dv->z_name, 0);
 748
 749           label = xasprintf (_("Z-score of %s"),var_to_string (dv->v));
 750           var_set_label (dst_var, label, false);
 751           free (label);
 752
 753           z = &t->z_scores[cnt++];
 754           z->src_var = dv->v;
 755           z->z_var = dst_var;
 756         }
 757     }
 758
 759   add_transformation (ds,
 760                       descriptives_trns_proc, descriptives_trns_free, t);
 761 }
 762 \f
 763 /* Statistical calculation. */
 764
 765 static bool listwise_missing (struct dsc_proc *dsc, const struct ccase *c);
 766
 767 /* Calculates and displays descriptive statistics for the cases
 768    in CF. */
 769 static void
 770 calc_descriptives (struct dsc_proc *dsc, struct casereader *group,
 771                    struct dataset *ds)
 772 {
 773   struct casereader *pass1, *pass2;
 774   casenumber count;
 775   struct ccase *c;
 776   size_t z_idx;
 777   size_t i;
 778
 779   c = casereader_peek (group, 0);
 780   if (c == NULL)
 781     {
 782       casereader_destroy (group);
 783       return;
 784     }
 785   output_split_file_values (ds, c);
 786   case_unref (c);
 787
 788   group = casereader_create_filter_weight (group, dataset_dict (ds),
 789                                            NULL, NULL);
 790
 791   pass1 = group;
 792   pass2 = dsc->max_moment <= MOMENT_MEAN ? NULL : casereader_clone (pass1);
 793
 794   for (i = 0; i < dsc->var_cnt; i++)
 795     {
 796       struct dsc_var *dv = &dsc->vars[i];
 797
 798       dv->valid = dv->missing = 0.0;
 799       if (dv->moments != NULL)
 800         moments_clear (dv->moments);
 801       dv->min = DBL_MAX;
 802       dv->max = -DBL_MAX;
 803     }
 804   dsc->missing_listwise = 0.;
 805   dsc->valid = 0.;
 806
 807   /* First pass to handle most of the work. */
 808   count = 0;
 809   for (; (c = casereader_read (pass1)) != NULL; case_unref (c))
 810     {
 811       double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
 812
 813       /* Check for missing values. */
 814       if (listwise_missing (dsc, c))
 815         {
 816           dsc->missing_listwise += weight;
 817           if (dsc->missing_type == DSC_LISTWISE)
 818             continue;
 819         }
 820       dsc->valid += weight;
 821
 822       for (i = 0; i < dsc->var_cnt; i++)
 823         {
 824           struct dsc_var *dv = &dsc->vars[i];
 825           double x = case_num (c, dv->v);
 826
 827           if (var_is_num_missing (dv->v, x, dsc->exclude))
 828             {
 829               dv->missing += weight;
 830               continue;
 831             }
 832
 833           if (dv->moments != NULL)
 834             moments_pass_one (dv->moments, x, weight);
 835
 836           if (x < dv->min)
 837             dv->min = x;
 838           if (x > dv->max)
 839             dv->max = x;
 840         }
 841
 842       count++;
 843     }
 844   if (!casereader_destroy (pass1))
 845     {
 846       casereader_destroy (pass2);
 847       return;
 848     }
 849
 850   /* Second pass for higher-order moments. */
 851   if (dsc->max_moment > MOMENT_MEAN)
 852     {
 853       for (; (c = casereader_read (pass2)) != NULL; case_unref (c))
 854         {
 855           double weight = dict_get_case_weight (dataset_dict (ds), c, NULL);
 856
 857           /* Check for missing values. */
 858           if (dsc->missing_type == DSC_LISTWISE && listwise_missing (dsc, c))
 859             continue;
 860
 861           for (i = 0; i < dsc->var_cnt; i++)
 862             {
 863               struct dsc_var *dv = &dsc->vars[i];
 864               double x = case_num (c, dv->v);
 865
 866               if (var_is_num_missing (dv->v, x, dsc->exclude))
 867                 continue;
 868
 869               if (dv->moments != NULL)
 870                 moments_pass_two (dv->moments, x, weight);
 871             }
 872         }
 873       if (!casereader_destroy (pass2))
 874         return;
 875     }
 876
 877   /* Calculate results. */
 878   if (dsc->z_writer)
 879     {
 880       c = case_create (casewriter_get_proto (dsc->z_writer));
 881       z_idx = 0;
 882       case_data_rw_idx (c, z_idx++)->f = count;
 883     }
 884   else
 885     c = NULL;
 886
 887   for (i = 0; i < dsc->var_cnt; i++)
 888     {
 889       struct dsc_var *dv = &dsc->vars[i];
 890       double W;
 891       int j;
 892
 893       for (j = 0; j < DSC_N_STATS; j++)
 894         dv->stats[j] = SYSMIS;
 895
 896       dv->valid = W = dsc->valid - dv->missing;
 897
 898       if (dv->moments != NULL)
 899         moments_calculate (dv->moments, NULL,
 900                            &dv->stats[DSC_MEAN], &dv->stats[DSC_VARIANCE],
 901                            &dv->stats[DSC_SKEWNESS], &dv->stats[DSC_KURTOSIS]);
 902       if (dsc->calc_stats & (1ul << DSC_SEMEAN)
 903           && dv->stats[DSC_VARIANCE] != SYSMIS && W > 0.)
 904         dv->stats[DSC_SEMEAN] = sqrt (dv->stats[DSC_VARIANCE]) / sqrt (W);
 905       if (dsc->calc_stats & (1ul << DSC_STDDEV)
 906           && dv->stats[DSC_VARIANCE] != SYSMIS)
 907         dv->stats[DSC_STDDEV] = sqrt (dv->stats[DSC_VARIANCE]);
 908       if (dsc->calc_stats & (1ul << DSC_SEKURT))
 909         if (dv->stats[DSC_KURTOSIS] != SYSMIS)
 910             dv->stats[DSC_SEKURT] = calc_sekurt (W);
 911       if (dsc->calc_stats & (1ul << DSC_SESKEW)
 912           && dv->stats[DSC_SKEWNESS] != SYSMIS)
 913         dv->stats[DSC_SESKEW] = calc_seskew (W);
 914       dv->stats[DSC_RANGE] = ((dv->min == DBL_MAX || dv->max == -DBL_MAX)
 915                               ? SYSMIS : dv->max - dv->min);
 916       dv->stats[DSC_MIN] = dv->min == DBL_MAX ? SYSMIS : dv->min;
 917       dv->stats[DSC_MAX] = dv->max == -DBL_MAX ? SYSMIS : dv->max;
 918       if (dsc->calc_stats & (1ul << DSC_SUM))
 919         dv->stats[DSC_SUM] = W * dv->stats[DSC_MEAN];
 920
 921       if (dv->z_name)
 922         {
 923           case_data_rw_idx (c, z_idx++)->f = dv->stats[DSC_MEAN];
 924           case_data_rw_idx (c, z_idx++)->f = dv->stats[DSC_STDDEV];
 925         }
 926     }
 927
 928   if (c != NULL)
 929     casewriter_write (dsc->z_writer, c);
 930
 931   /* Output results. */
 932   display (dsc);
 933 }
 934
 935 /* Returns true if any of the descriptives variables in DSC's
 936    variable list have missing values in case C, false otherwise. */
 937 static bool
 938 listwise_missing (struct dsc_proc *dsc, const struct ccase *c)
 939 {
 940   size_t i;
 941
 942   for (i = 0; i < dsc->var_cnt; i++)
 943     {
 944       struct dsc_var *dv = &dsc->vars[i];
 945       double x = case_num (c, dv->v);
 946
 947       if (var_is_num_missing (dv->v, x, dsc->exclude))
 948         return true;
 949     }
 950   return false;
 951 }
 952 \f
 953 /* Statistical display. */
 954
 955 static algo_compare_func descriptives_compare_dsc_vars;
 956
 957 /* Displays a table of descriptive statistics for DSC. */
 958 static void
 959 display (struct dsc_proc *dsc)
 960 {
 961   size_t i;
 962   int nc;
 963   struct tab_table *t;
 964
 965   nc = 1 + (dsc->format == DSC_SERIAL ? 2 : 1);
 966   for (i = 0; i < DSC_N_STATS; i++)
 967     if (dsc->show_stats & (1ul << i))
 968       nc++;
 969
 970   if (dsc->sort_by_stat != DSC_NONE)
 971     sort (dsc->vars, dsc->var_cnt, sizeof *dsc->vars,
 972           descriptives_compare_dsc_vars, dsc);
 973
 974   t = tab_create (nc, dsc->var_cnt + 1);
 975   tab_headers (t, 1, 0, 1, 0);
 976   tab_box (t, TAL_1, TAL_1, -1, -1, 0, 0, nc - 1, dsc->var_cnt);
 977   tab_box (t, -1, -1, -1, TAL_1, 1, 0, nc - 1, dsc->var_cnt);
 978   tab_hline (t, TAL_2, 0, nc - 1, 1);
 979   tab_vline (t, TAL_2, 1, 0, dsc->var_cnt);
 980
 981   nc = 0;
 982   tab_text (t, nc++, 0, TAB_LEFT | TAT_TITLE, _("Variable"));
 983   if (dsc->format == DSC_SERIAL)
 984     {
 985       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Valid N"));
 986       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Missing N"));
 987     }
 988   else
 989     tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, "N");
 990
 991   for (i = 0; i < DSC_N_STATS; i++)
 992     if (dsc->show_stats & (1ul << i))
 993       {
 994         const char *title = gettext (dsc_info[i].name);
 995         tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, title);
 996       }
 997
 998   for (i = 0; i < dsc->var_cnt; i++)
 999     {
1000       struct dsc_var *dv = &dsc->vars[i];
1001       size_t j;
1002
1003       nc = 0;
1004       tab_text (t, nc++, i + 1, TAB_LEFT, var_to_string (dv->v));
1005       tab_text_format (t, nc++, i + 1, 0, "%g", dv->valid);
1006       if (dsc->format == DSC_SERIAL)
1007         tab_text_format (t, nc++, i + 1, 0, "%g", dv->missing);
1008
1009       for (j = 0; j < DSC_N_STATS; j++)
1010         if (dsc->show_stats & (1ul << j))
1011           tab_double (t, nc++, i + 1, TAB_NONE, dv->stats[j], NULL);
1012     }
1013
1014   tab_title (t, _("Valid cases = %g; cases with missing value(s) = %g."),
1015              dsc->valid, dsc->missing_listwise);
1016
1017   tab_submit (t);
1018 }
1019
1020 /* Compares `struct dsc_var's A and B according to the ordering
1021    specified by CMD. */
1022 static int
1023 descriptives_compare_dsc_vars (const void *a_, const void *b_, const void *dsc_)
1024 {
1025   const struct dsc_var *a = a_;
1026   const struct dsc_var *b = b_;
1027   const struct dsc_proc *dsc = dsc_;
1028
1029   int result;
1030
1031   if (dsc->sort_by_stat == DSC_NAME)
1032     result = utf8_strcasecmp (var_get_name (a->v), var_get_name (b->v));
1033   else
1034     {
1035       double as = a->stats[dsc->sort_by_stat];
1036       double bs = b->stats[dsc->sort_by_stat];
1037
1038       result = as < bs ? -1 : as > bs;
1039     }
1040
1041   if (!dsc->sort_ascending)
1042     result = -result;
1043
1044   return result;
1045 }