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