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