src/language/stats/descriptives.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
   3    Written by Ben Pfaff <blp@gnu.org>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful, but
  11    WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    General Public License for more details.
  14
  15    You should have received a copy of the GNU General Public License
  16    along with this program; if not, write to the Free Software
  17    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  18    02110-1301, USA. */
  19
  20 /* FIXME: Many possible optimizations. */
  21
  22 #include <config.h>
  23 #include <libpspp/message.h>
  24 #include <limits.h>
  25 #include <math.h>
  26 #include <stdlib.h>
  27 #include <libpspp/array.h>
  28 #include <libpspp/alloc.h>
  29 #include <data/case.h>
  30 #include <data/casefile.h>
  31 #include <language/command.h>
  32 #include <libpspp/compiler.h>
  33 #include <data/dictionary.h>
  34 #include <language/lexer/lexer.h>
  35 #include <libpspp/message.h>
  36 #include <libpspp/magic.h>
  37 #include <math/moments.h>
  38 #include <output/manager.h>
  39 #include <output/table.h>
  40 #include <data/variable.h>
  41 #include <procedure.h>
  42
  43 #include "gettext.h"
  44 #define _(msgid) gettext (msgid)
  45 #define N_(msgid) msgid
  46
  47 /* DESCRIPTIVES private data. */
  48
  49 struct dsc_proc;
  50
  51 /* Handling of missing values. */
  52 enum dsc_missing_type
  53   {
  54     DSC_VARIABLE,       /* Handle missing values on a per-variable basis. */
  55     DSC_LISTWISE        /* Discard entire case if any variable is missing. */
  56   };
  57
  58 /* Describes properties of a distribution for the purpose of
  59    calculating a Z-score. */
  60 struct dsc_z_score
  61   {
  62     int src_idx;                /* Source index into case data. */
  63     int dst_idx;                /* Destination index into case data. */
  64     double mean;                /* Distribution mean. */
  65     double std_dev;             /* Distribution standard deviation. */
  66     struct variable *v;         /* Variable on which z-score is based. */
  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     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     int include_user_missing;   /* Nonzero to include user-missing values. */
  78   };
  79
  80 /* Statistics.  Used as bit indexes, so must be 32 or fewer. */
  81 enum dsc_statistic
  82   {
  83     DSC_MEAN = 0, DSC_SEMEAN, DSC_STDDEV, DSC_VARIANCE, DSC_KURTOSIS,
  84     DSC_SEKURT, DSC_SKEWNESS, DSC_SESKEW, DSC_RANGE, DSC_MIN,
  85     DSC_MAX, DSC_SUM, DSC_N_STATS,
  86
  87     /* Only valid as sort criteria. */
  88     DSC_NAME = -2,              /* Sort by name. */
  89     DSC_NONE = -1               /* Unsorted. */
  90   };
  91
  92 /* Describes one statistic. */
  93 struct dsc_statistic_info
  94   {
  95     const char *identifier;     /* Identifier. */
  96     const char *name;           /* Full name. */
  97     enum moment moment;         /* Highest moment needed to calculate. */
  98   };
  99
 100 /* Table of statistics, indexed by DSC_*. */
 101 static const struct dsc_statistic_info dsc_info[DSC_N_STATS] =
 102   {
 103     {"MEAN", N_("Mean"), MOMENT_MEAN},
 104     {"SEMEAN", N_("S E Mean"), MOMENT_VARIANCE},
 105     {"STDDEV", N_("Std Dev"), MOMENT_VARIANCE},
 106     {"VARIANCE", N_("Variance"), MOMENT_VARIANCE},
 107     {"KURTOSIS", N_("Kurtosis"), MOMENT_KURTOSIS},
 108     {"SEKURTOSIS", N_("S E Kurt"), MOMENT_NONE},
 109     {"SKEWNESS", N_("Skewness"), MOMENT_SKEWNESS},
 110     {"SESKEWNESS", N_("S E Skew"), MOMENT_NONE},
 111     {"RANGE", N_("Range"), MOMENT_NONE},
 112     {"MINIMUM", N_("Minimum"), MOMENT_NONE},
 113     {"MAXIMUM", N_("Maximum"), MOMENT_NONE},
 114     {"SUM", N_("Sum"), MOMENT_MEAN},
 115   };
 116
 117 /* Statistics calculated by default if none are explicitly
 118    requested. */
 119 #define DEFAULT_STATS                                                   \
 120         ((1ul << DSC_MEAN) | (1ul << DSC_STDDEV) | (1ul << DSC_MIN)     \
 121          | (1ul << DSC_MAX))
 122
 123 /* A variable specified on DESCRIPTIVES. */
 124 struct dsc_var
 125   {
 126     struct variable *v;         /* Variable to calculate on. */
 127     char z_name[LONG_NAME_LEN + 1]; /* Name for z-score variable. */
 128     double valid, missing;      /* Valid, missing counts. */
 129     struct moments *moments;    /* Moments. */
 130     double min, max;            /* Maximum and mimimum values. */
 131     double stats[DSC_N_STATS];  /* All the stats' values. */
 132   };
 133
 134 /* Output format. */
 135 enum dsc_format
 136   {
 137     DSC_LINE,           /* Abbreviated format. */
 138     DSC_SERIAL          /* Long format. */
 139   };
 140
 141 /* A DESCRIPTIVES procedure. */
 142 struct dsc_proc
 143   {
 144     /* Per-variable info. */
 145     struct dsc_var *vars;       /* Variables. */
 146     size_t var_cnt;             /* Number of variables. */
 147
 148     /* User options. */
 149     enum dsc_missing_type missing_type; /* Treatment of missing values. */
 150     int include_user_missing;   /* Nonzero to include user-missing values. */
 151     int show_var_labels;        /* Nonzero to show variable labels. */
 152     int show_index;             /* Nonzero to show variable index. */
 153     enum dsc_format format;     /* Output format. */
 154
 155     /* Accumulated results. */
 156     double missing_listwise;    /* Sum of weights of cases missing listwise. */
 157     double valid;               /* Sum of weights of valid cases. */
 158     int bad_warn;               /* Warn if bad weight found. */
 159     enum dsc_statistic sort_by_stat; /* Statistic to sort by; -1: name. */
 160     int sort_ascending;         /* !0: ascending order; 0: descending. */
 161     unsigned long show_stats;   /* Statistics to display. */
 162     unsigned long calc_stats;   /* Statistics to calculate. */
 163     enum moment max_moment;     /* Highest moment needed for stats. */
 164   };
 165
 166 /* Parsing. */
 167 static enum dsc_statistic match_statistic (void);
 168 static void free_dsc_proc (struct dsc_proc *);
 169
 170 /* Z-score functions. */
 171 static int try_name (struct dsc_proc *dsc, char *name);
 172 static int generate_z_varname (struct dsc_proc *dsc, char *z_name,
 173                                const char *name, size_t *z_cnt);
 174 static void dump_z_table (struct dsc_proc *);
 175 static void setup_z_trns (struct dsc_proc *);
 176
 177 /* Procedure execution functions. */
 178 static bool calc_descriptives (const struct casefile *, void *dsc_);
 179 static void display (struct dsc_proc *dsc);
 180 \f
 181 /* Parser and outline. */
 182
 183 /* Handles DESCRIPTIVES. */
 184 int
 185 cmd_descriptives (void)
 186 {
 187   struct dsc_proc *dsc;
 188   struct variable **vars = NULL;
 189   size_t var_cnt = 0;
 190   int save_z_scores = 0;
 191   size_t z_cnt = 0;
 192   size_t i;
 193   bool ok;
 194
 195   /* Create and initialize dsc. */
 196   dsc = xmalloc (sizeof *dsc);
 197   dsc->vars = NULL;
 198   dsc->var_cnt = 0;
 199   dsc->missing_type = DSC_VARIABLE;
 200   dsc->include_user_missing = 0;
 201   dsc->show_var_labels = 1;
 202   dsc->show_index = 0;
 203   dsc->format = DSC_LINE;
 204   dsc->missing_listwise = 0.;
 205   dsc->valid = 0.;
 206   dsc->bad_warn = 1;
 207   dsc->sort_by_stat = DSC_NONE;
 208   dsc->sort_ascending = 1;
 209   dsc->show_stats = dsc->calc_stats = DEFAULT_STATS;
 210
 211   /* Parse DESCRIPTIVES. */
 212   while (token != '.')
 213     {
 214       if (lex_match_id ("MISSING"))
 215         {
 216           lex_match ('=');
 217           while (token != '.' && token != '/')
 218             {
 219               if (lex_match_id ("VARIABLE"))
 220                 dsc->missing_type = DSC_VARIABLE;
 221               else if (lex_match_id ("LISTWISE"))
 222                 dsc->missing_type = DSC_LISTWISE;
 223               else if (lex_match_id ("INCLUDE"))
 224                 dsc->include_user_missing = 1;
 225               else
 226                 {
 227                   lex_error (NULL);
 228                   goto error;
 229                 }
 230               lex_match (',');
 231             }
 232         }
 233       else if (lex_match_id ("SAVE"))
 234         save_z_scores = 1;
 235       else if (lex_match_id ("FORMAT"))
 236         {
 237           lex_match ('=');
 238           while (token != '.' && token != '/')
 239             {
 240               if (lex_match_id ("LABELS"))
 241                 dsc->show_var_labels = 1;
 242               else if (lex_match_id ("NOLABELS"))
 243                 dsc->show_var_labels = 0;
 244               else if (lex_match_id ("INDEX"))
 245                 dsc->show_index = 1;
 246               else if (lex_match_id ("NOINDEX"))
 247                 dsc->show_index = 0;
 248               else if (lex_match_id ("LINE"))
 249                 dsc->format = DSC_LINE;
 250               else if (lex_match_id ("SERIAL"))
 251                 dsc->format = DSC_SERIAL;
 252               else
 253                 {
 254                   lex_error (NULL);
 255                   goto error;
 256                 }
 257               lex_match (',');
 258             }
 259         }
 260       else if (lex_match_id ("STATISTICS"))
 261         {
 262           lex_match ('=');
 263           dsc->show_stats = 0;
 264           while (token != '.' && token != '/')
 265             {
 266               if (lex_match (T_ALL))
 267                 dsc->show_stats |= (1ul << DSC_N_STATS) - 1;
 268               else if (lex_match_id ("DEFAULT"))
 269                 dsc->show_stats |= DEFAULT_STATS;
 270               else
 271                 dsc->show_stats |= 1ul << (match_statistic ());
 272               lex_match (',');
 273             }
 274           if (dsc->show_stats == 0)
 275             dsc->show_stats = DEFAULT_STATS;
 276         }
 277       else if (lex_match_id ("SORT"))
 278         {
 279           lex_match ('=');
 280           if (lex_match_id ("NAME"))
 281             dsc->sort_by_stat = DSC_NAME;
 282           else
 283             {
 284               dsc->sort_by_stat = match_statistic ();
 285               if (dsc->sort_by_stat == DSC_NONE )
 286                 dsc->sort_by_stat = DSC_MEAN;
 287             }
 288           if (lex_match ('('))
 289             {
 290               if (lex_match_id ("A"))
 291                 dsc->sort_ascending = 1;
 292               else if (lex_match_id ("D"))
 293                 dsc->sort_ascending = 0;
 294               else
 295                 lex_error (NULL);
 296               lex_force_match (')');
 297             }
 298         }
 299       else if (var_cnt == 0)
 300         {
 301           if (lex_look_ahead () == '=')
 302             {
 303               lex_match_id ("VARIABLES");
 304               lex_match ('=');
 305             }
 306
 307           while (token != '.' && token != '/')
 308             {
 309               int i;
 310
 311               if (!parse_variables (default_dict, &vars, &var_cnt,
 312                                     PV_APPEND | PV_NO_DUPLICATE | PV_NUMERIC))
 313                 goto error;
 314
 315               dsc->vars = xnrealloc (dsc->vars, var_cnt, sizeof *dsc->vars);
 316               for (i = dsc->var_cnt; i < var_cnt; i++)
 317                 {
 318                   struct dsc_var *dv = &dsc->vars[i];
 319                   dv->v = vars[i];
 320                   dv->z_name[0] = '\0';
 321                   dv->moments = NULL;
 322                 }
 323               dsc->var_cnt = var_cnt;
 324
 325               if (lex_match ('('))
 326                 {
 327                   if (token != T_ID)
 328                     {
 329                       lex_error (NULL);
 330                       goto error;
 331                     }
 332                   if (try_name (dsc, tokid))
 333                     {
 334                       strcpy (dsc->vars[dsc->var_cnt - 1].z_name, tokid);
 335                       z_cnt++;
 336                     }
 337                   else
 338                     msg (SE, _("Z-score variable name %s would be"
 339                                " a duplicate variable name."), tokid);
 340                   lex_get ();
 341                   if (!lex_force_match (')'))
 342                     goto error;
 343                 }
 344             }
 345         }
 346       else
 347         {
 348           lex_error (NULL);
 349           goto error;
 350         }
 351
 352       lex_match ('/');
 353     }
 354   if (var_cnt == 0)
 355     {
 356       msg (SE, _("No variables specified."));
 357       goto error;
 358     }
 359
 360   /* Construct z-score varnames, show translation table. */
 361   if (z_cnt || save_z_scores)
 362     {
 363       if (save_z_scores)
 364         {
 365           size_t gen_cnt = 0;
 366
 367           for (i = 0; i < dsc->var_cnt; i++)
 368             if (dsc->vars[i].z_name[0] == 0)
 369               {
 370                 if (!generate_z_varname (dsc, dsc->vars[i].z_name,
 371                                          dsc->vars[i].v->name, &gen_cnt))
 372                   goto error;
 373                 z_cnt++;
 374               }
 375         }
 376       dump_z_table (dsc);
 377     }
 378
 379   /* Figure out statistics to display. */
 380   if (dsc->show_stats & (1ul << DSC_SKEWNESS))
 381     dsc->show_stats |= 1ul << DSC_SESKEW;
 382   if (dsc->show_stats & (1ul << DSC_KURTOSIS))
 383     dsc->show_stats |= 1ul << DSC_SEKURT;
 384
 385   /* Figure out which statistics to calculate. */
 386   dsc->calc_stats = dsc->show_stats;
 387   if (z_cnt > 0)
 388     dsc->calc_stats |= (1ul << DSC_MEAN) | (1ul << DSC_STDDEV);
 389   if (dsc->sort_by_stat >= 0)
 390     dsc->calc_stats |= 1ul << dsc->sort_by_stat;
 391   if (dsc->show_stats & (1ul << DSC_SESKEW))
 392     dsc->calc_stats |= 1ul << DSC_SKEWNESS;
 393   if (dsc->show_stats & (1ul << DSC_SEKURT))
 394     dsc->calc_stats |= 1ul << DSC_KURTOSIS;
 395
 396   /* Figure out maximum moment needed and allocate moments for
 397      the variables. */
 398   dsc->max_moment = MOMENT_NONE;
 399   for (i = 0; i < DSC_N_STATS; i++)
 400     if (dsc->calc_stats & (1ul << i) && dsc_info[i].moment > dsc->max_moment)
 401       dsc->max_moment = dsc_info[i].moment;
 402   if (dsc->max_moment != MOMENT_NONE)
 403     for (i = 0; i < dsc->var_cnt; i++)
 404       dsc->vars[i].moments = moments_create (dsc->max_moment);
 405
 406   /* Data pass. */
 407   ok = multipass_procedure_with_splits (calc_descriptives, dsc);
 408
 409   /* Z-scoring! */
 410   if (ok && z_cnt)
 411     setup_z_trns (dsc);
 412
 413   /* Done. */
 414   free (vars);
 415   free_dsc_proc (dsc);
 416   return ok ? CMD_SUCCESS : CMD_CASCADING_FAILURE;
 417
 418  error:
 419   free (vars);
 420   free_dsc_proc (dsc);
 421   return CMD_FAILURE;
 422 }
 423
 424 /* Returns the statistic named by the current token and skips past the token.
 425    Returns DSC_NONE if no statistic is given (e.g., subcommand with no
 426    specifiers). Emits an error if the current token ID does not name a
 427    statistic. */
 428 static enum dsc_statistic
 429 match_statistic (void)
 430 {
 431   if (token == T_ID)
 432     {
 433       enum dsc_statistic stat;
 434
 435       for (stat = 0; stat < DSC_N_STATS; stat++)
 436         if (lex_match_id (dsc_info[stat].identifier))
 437           return stat;
 438
 439       lex_get();
 440       lex_error (_("expecting statistic name: reverting to default"));
 441     }
 442
 443   return DSC_NONE;
 444 }
 445
 446 /* Frees DSC. */
 447 static void
 448 free_dsc_proc (struct dsc_proc *dsc)
 449 {
 450   size_t i;
 451
 452   if (dsc == NULL)
 453     return;
 454
 455   for (i = 0; i < dsc->var_cnt; i++)
 456     moments_destroy (dsc->vars[i].moments);
 457   free (dsc->vars);
 458   free (dsc);
 459 }
 460 \f
 461 /* Z scores. */
 462
 463 /* Returns 0 if NAME is a duplicate of any existing variable name or
 464    of any previously-declared z-var name; otherwise returns 1. */
 465 static int
 466 try_name (struct dsc_proc *dsc, char *name)
 467 {
 468   size_t i;
 469
 470   if (dict_lookup_var (default_dict, name) != NULL)
 471     return 0;
 472   for (i = 0; i < dsc->var_cnt; i++)
 473     if (!strcasecmp (dsc->vars[i].z_name, name))
 474       return 0;
 475   return 1;
 476 }
 477
 478 /* Generates a name for a Z-score variable based on a variable
 479    named VAR_NAME, given that *Z_CNT generated variable names are
 480    known to already exist.  If successful, returns nonzero and
 481    copies the new name into Z_NAME.  On failure, returns zero. */
 482 static int
 483 generate_z_varname (struct dsc_proc *dsc, char *z_name,
 484                     const char *var_name, size_t *z_cnt)
 485 {
 486   char name[LONG_NAME_LEN + 1];
 487
 488   /* Try a name based on the original variable name. */
 489   name[0] = 'Z';
 490   str_copy_trunc (name + 1, sizeof name - 1, var_name);
 491   if (try_name (dsc, name))
 492     {
 493       strcpy (z_name, name);
 494       return 1;
 495     }
 496
 497   /* Generate a synthetic name. */
 498   for (;;)
 499     {
 500       (*z_cnt)++;
 501
 502       if (*z_cnt <= 99)
 503         sprintf (name, "ZSC%03d", *z_cnt);
 504       else if (*z_cnt <= 108)
 505         sprintf (name, "STDZ%02d", *z_cnt - 99);
 506       else if (*z_cnt <= 117)
 507         sprintf (name, "ZZZZ%02d", *z_cnt - 108);
 508       else if (*z_cnt <= 126)
 509         sprintf (name, "ZQZQ%02d", *z_cnt - 117);
 510       else
 511         {
 512           msg (SE, _("Ran out of generic names for Z-score variables.  "
 513                      "There are only 126 generic names: ZSC001-ZSC0999, "
 514                      "STDZ01-STDZ09, ZZZZ01-ZZZZ09, ZQZQ01-ZQZQ09."));
 515           return 0;
 516         }
 517
 518       if (try_name (dsc, name))
 519         {
 520           strcpy (z_name, name);
 521           return 1;
 522         }
 523     }
 524 }
 525
 526 /* Outputs a table describing the mapping between source
 527    variables and Z-score variables. */
 528 static void
 529 dump_z_table (struct dsc_proc *dsc)
 530 {
 531   size_t cnt = 0;
 532   struct tab_table *t;
 533
 534   {
 535     size_t i;
 536
 537     for (i = 0; i < dsc->var_cnt; i++)
 538       if (dsc->vars[i].z_name[0] != '\0')
 539         cnt++;
 540   }
 541
 542   t = tab_create (2, cnt + 1, 0);
 543   tab_title (t, 0, _("Mapping of variables to corresponding Z-scores."));
 544   tab_columns (t, SOM_COL_DOWN, 1);
 545   tab_headers (t, 0, 0, 1, 0);
 546   tab_box (t, TAL_1, TAL_1, TAL_0, TAL_1, 0, 0, 1, cnt);
 547   tab_hline (t, TAL_2, 0, 1, 1);
 548   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
 549   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Target"));
 550   tab_dim (t, tab_natural_dimensions);
 551
 552   {
 553     size_t i, y;
 554
 555     for (i = 0, y = 1; i < dsc->var_cnt; i++)
 556       if (dsc->vars[i].z_name[0] != '\0')
 557         {
 558           tab_text (t, 0, y, TAB_LEFT, dsc->vars[i].v->name);
 559           tab_text (t, 1, y++, TAB_LEFT, dsc->vars[i].z_name);
 560         }
 561   }
 562
 563   tab_submit (t);
 564 }
 565
 566 /* Transformation function to calculate Z-scores. Will return SYSMIS if any of
 567    the following are true: 1) mean or standard deviation is SYSMIS 2) score is
 568    SYSMIS 3) score is user missing and they were not included in the original
 569    analyis. 4) any of the variables in the original analysis were missing
 570    (either system or user-missing values that weren't included).
 571 */
 572 static int
 573 descriptives_trns_proc (void *trns_, struct ccase * c,
 574                         int case_idx UNUSED)
 575 {
 576   struct dsc_trns *t = trns_;
 577   struct dsc_z_score *z;
 578   struct variable **vars;
 579   int all_sysmis = 0;
 580
 581   if (t->missing_type == DSC_LISTWISE)
 582     {
 583       assert(t->vars);
 584       for (vars = t->vars; vars < t->vars + t->var_cnt; vars++)
 585         {
 586           double score = case_num (c, (*vars)->fv);
 587           if ( score == SYSMIS
 588                || (!t->include_user_missing
 589                    && mv_is_num_user_missing (&(*vars)->miss, score)))
 590             {
 591               all_sysmis = 1;
 592               break;
 593             }
 594         }
 595     }
 596
 597   for (z = t->z_scores; z < t->z_scores + t->z_score_cnt; z++)
 598     {
 599       double input = case_num (c, z->src_idx);
 600       double *output = &case_data_rw (c, z->dst_idx)->f;
 601
 602       if (z->mean == SYSMIS || z->std_dev == SYSMIS
 603           || all_sysmis || input == SYSMIS
 604           || (!t->include_user_missing
 605               && mv_is_num_user_missing (&z->v->miss, input)))
 606         *output = SYSMIS;
 607       else
 608         *output = (input - z->mean) / z->std_dev;
 609     }
 610   return TRNS_CONTINUE;
 611 }
 612
 613 /* Frees a descriptives_trns struct. */
 614 static bool
 615 descriptives_trns_free (void *trns_)
 616 {
 617   struct dsc_trns *t = trns_;
 618
 619   free (t->z_scores);
 620   assert((t->missing_type != DSC_LISTWISE) ^ (t->vars != NULL));
 621   free (t->vars);
 622   return true;
 623 }
 624
 625 /* Sets up a transformation to calculate Z scores. */
 626 static void
 627 setup_z_trns (struct dsc_proc *dsc)
 628 {
 629   struct dsc_trns *t;
 630   size_t cnt, i;
 631
 632   for (cnt = i = 0; i < dsc->var_cnt; i++)
 633     if (dsc->vars[i].z_name[0] != '\0')
 634       cnt++;
 635
 636   t = xmalloc (sizeof *t);
 637   t->z_scores = xnmalloc (cnt, sizeof *t->z_scores);
 638   t->z_score_cnt = cnt;
 639   t->missing_type = dsc->missing_type;
 640   t->include_user_missing = dsc->include_user_missing;
 641   if ( t->missing_type == DSC_LISTWISE )
 642     {
 643       t->var_cnt = dsc->var_cnt;
 644       t->vars = xnmalloc (t->var_cnt, sizeof *t->vars);
 645       for (i = 0; i < t->var_cnt; i++)
 646         t->vars[i] = dsc->vars[i].v;
 647     }
 648   else
 649     {
 650       t->var_cnt = 0;
 651       t->vars = NULL;
 652     }
 653
 654   for (cnt = i = 0; i < dsc->var_cnt; i++)
 655     {
 656       struct dsc_var *dv = &dsc->vars[i];
 657       if (dv->z_name[0] != '\0')
 658         {
 659           struct dsc_z_score *z;
 660           char *cp;
 661           struct variable *dst_var;
 662
 663           dst_var = dict_create_var_assert (default_dict, dv->z_name, 0);
 664           dst_var->init = 0;
 665           if (dv->v->label)
 666             {
 667               dst_var->label = xmalloc (strlen (dv->v->label) + 12);
 668               cp = stpcpy (dst_var->label, _("Z-score of "));
 669               strcpy (cp, dv->v->label);
 670             }
 671           else
 672             {
 673               dst_var->label = xmalloc (strlen (dv->v->name) + 12);
 674               cp = stpcpy (dst_var->label, _("Z-score of "));
 675               strcpy (cp, dv->v->name);
 676             }
 677
 678           z = &t->z_scores[cnt++];
 679           z->src_idx = dv->v->fv;
 680           z->dst_idx = dst_var->fv;
 681           z->mean = dv->stats[DSC_MEAN];
 682           z->std_dev = dv->stats[DSC_STDDEV];
 683           z->v = dv->v;
 684         }
 685     }
 686
 687   add_transformation (descriptives_trns_proc, descriptives_trns_free, t);
 688 }
 689 \f
 690 /* Statistical calculation. */
 691
 692 static int listwise_missing (struct dsc_proc *dsc, const struct ccase *c);
 693
 694 /* Calculates and displays descriptive statistics for the cases
 695    in CF. */
 696 static bool
 697 calc_descriptives (const struct casefile *cf, void *dsc_)
 698 {
 699   struct dsc_proc *dsc = dsc_;
 700   struct casereader *reader;
 701   struct ccase c;
 702   size_t i;
 703
 704   for (i = 0; i < dsc->var_cnt; i++)
 705     {
 706       struct dsc_var *dv = &dsc->vars[i];
 707
 708       dv->valid = dv->missing = 0.0;
 709       if (dv->moments != NULL)
 710         moments_clear (dv->moments);
 711       dv->min = DBL_MAX;
 712       dv->max = -DBL_MAX;
 713     }
 714   dsc->missing_listwise = 0.;
 715   dsc->valid = 0.;
 716
 717   /* First pass to handle most of the work. */
 718   for (reader = casefile_get_reader (cf);
 719        casereader_read (reader, &c);
 720        case_destroy (&c))
 721     {
 722       double weight = dict_get_case_weight (default_dict, &c, &dsc->bad_warn);
 723       if (weight <= 0.0)
 724         continue;
 725
 726       /* Check for missing values. */
 727       if (listwise_missing (dsc, &c))
 728         {
 729           dsc->missing_listwise += weight;
 730           if (dsc->missing_type == DSC_LISTWISE)
 731             continue;
 732         }
 733       dsc->valid += weight;
 734
 735       for (i = 0; i < dsc->var_cnt; i++)
 736         {
 737           struct dsc_var *dv = &dsc->vars[i];
 738           double x = case_num (&c, dv->v->fv);
 739
 740           if (dsc->missing_type != DSC_LISTWISE
 741               && (x == SYSMIS
 742                   || (!dsc->include_user_missing
 743                       && mv_is_num_user_missing (&dv->v->miss, x))))
 744             {
 745               dv->missing += weight;
 746               continue;
 747             }
 748
 749           if (dv->moments != NULL)
 750             moments_pass_one (dv->moments, x, weight);
 751
 752           if (x < dv->min)
 753             dv->min = x;
 754           if (x > dv->max)
 755             dv->max = x;
 756         }
 757     }
 758   casereader_destroy (reader);
 759
 760   /* Second pass for higher-order moments. */
 761   if (dsc->max_moment > MOMENT_MEAN)
 762     {
 763       for (reader = casefile_get_reader (cf);
 764            casereader_read (reader, &c);
 765            case_destroy (&c))
 766         {
 767           double weight = dict_get_case_weight (default_dict, &c,
 768                                                 &dsc->bad_warn);
 769           if (weight <= 0.0)
 770             continue;
 771
 772           /* Check for missing values. */
 773           if (listwise_missing (dsc, &c)
 774               && dsc->missing_type == DSC_LISTWISE)
 775             continue;
 776
 777           for (i = 0; i < dsc->var_cnt; i++)
 778             {
 779               struct dsc_var *dv = &dsc->vars[i];
 780               double x = case_num (&c, dv->v->fv);
 781
 782               if (dsc->missing_type != DSC_LISTWISE
 783                   && (x == SYSMIS
 784                       || (!dsc->include_user_missing
 785                           && mv_is_num_user_missing (&dv->v->miss, x))))
 786                 continue;
 787
 788               if (dv->moments != NULL)
 789                 moments_pass_two (dv->moments, x, weight);
 790             }
 791         }
 792       casereader_destroy (reader);
 793     }
 794
 795   /* Calculate results. */
 796   for (i = 0; i < dsc->var_cnt; i++)
 797     {
 798       struct dsc_var *dv = &dsc->vars[i];
 799       double W;
 800       int j;
 801
 802       for (j = 0; j < DSC_N_STATS; j++)
 803         dv->stats[j] = SYSMIS;
 804
 805       dv->valid = W = dsc->valid - dv->missing;
 806
 807       if (dv->moments != NULL)
 808         moments_calculate (dv->moments, NULL,
 809                            &dv->stats[DSC_MEAN], &dv->stats[DSC_VARIANCE],
 810                            &dv->stats[DSC_SKEWNESS], &dv->stats[DSC_KURTOSIS]);
 811       if (dsc->calc_stats & (1ul << DSC_SEMEAN)
 812           && dv->stats[DSC_VARIANCE] != SYSMIS && W > 0.)
 813         dv->stats[DSC_SEMEAN] = sqrt (dv->stats[DSC_VARIANCE]) / sqrt (W);
 814       if (dsc->calc_stats & (1ul << DSC_STDDEV)
 815           && dv->stats[DSC_VARIANCE] != SYSMIS)
 816         dv->stats[DSC_STDDEV] = sqrt (dv->stats[DSC_VARIANCE]);
 817       if (dsc->calc_stats & (1ul << DSC_SEKURT))
 818         if (dv->stats[DSC_KURTOSIS] != SYSMIS)
 819             dv->stats[DSC_SEKURT] = calc_sekurt (W);
 820       if (dsc->calc_stats & (1ul << DSC_SESKEW)
 821           && dv->stats[DSC_SKEWNESS] != SYSMIS)
 822         dv->stats[DSC_SESKEW] = calc_seskew (W);
 823       dv->stats[DSC_RANGE] = ((dv->min == DBL_MAX || dv->max == -DBL_MAX)
 824                               ? SYSMIS : dv->max - dv->min);
 825       dv->stats[DSC_MIN] = dv->min == DBL_MAX ? SYSMIS : dv->min;
 826       dv->stats[DSC_MAX] = dv->max == -DBL_MAX ? SYSMIS : dv->max;
 827       if (dsc->calc_stats & (1ul << DSC_SUM))
 828         dv->stats[DSC_SUM] = W * dv->stats[DSC_MEAN];
 829     }
 830
 831   /* Output results. */
 832   display (dsc);
 833
 834   return true;
 835 }
 836
 837 /* Returns nonzero if any of the descriptives variables in DSC's
 838    variable list have missing values in case C, zero otherwise. */
 839 static int
 840 listwise_missing (struct dsc_proc *dsc, const struct ccase *c)
 841 {
 842   size_t i;
 843
 844   for (i = 0; i < dsc->var_cnt; i++)
 845     {
 846       struct dsc_var *dv = &dsc->vars[i];
 847       double x = case_num (c, dv->v->fv);
 848
 849       if (x == SYSMIS
 850           || (!dsc->include_user_missing
 851               && mv_is_num_user_missing (&dv->v->miss, x)))
 852         return 1;
 853     }
 854   return 0;
 855 }
 856 \f
 857 /* Statistical display. */
 858
 859 static algo_compare_func descriptives_compare_dsc_vars;
 860
 861 /* Displays a table of descriptive statistics for DSC. */
 862 static void
 863 display (struct dsc_proc *dsc)
 864 {
 865   size_t i;
 866   int nc;
 867   struct tab_table *t;
 868
 869   nc = 1 + (dsc->format == DSC_SERIAL ? 2 : 1);
 870   for (i = 0; i < DSC_N_STATS; i++)
 871     if (dsc->show_stats & (1ul << i))
 872       nc++;
 873
 874   if (dsc->sort_by_stat != DSC_NONE)
 875     sort (dsc->vars, dsc->var_cnt, sizeof *dsc->vars,
 876           descriptives_compare_dsc_vars, dsc);
 877
 878   t = tab_create (nc, dsc->var_cnt + 1, 0);
 879   tab_headers (t, 1, 0, 1, 0);
 880   tab_box (t, TAL_1, TAL_1, -1, -1, 0, 0, nc - 1, dsc->var_cnt);
 881   tab_box (t, -1, -1, -1, TAL_1, 1, 0, nc - 1, dsc->var_cnt);
 882   tab_hline (t, TAL_2, 0, nc - 1, 1);
 883   tab_vline (t, TAL_2, 1, 0, dsc->var_cnt);
 884   tab_dim (t, tab_natural_dimensions);
 885
 886   nc = 0;
 887   tab_text (t, nc++, 0, TAB_LEFT | TAT_TITLE, _("Variable"));
 888   if (dsc->format == DSC_SERIAL)
 889     {
 890       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Valid N"));
 891       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Missing N"));
 892     }
 893   else
 894     tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, "N");
 895
 896   for (i = 0; i < DSC_N_STATS; i++)
 897     if (dsc->show_stats & (1ul << i))
 898       {
 899         const char *title = gettext (dsc_info[i].name);
 900         tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, title);
 901       }
 902
 903   for (i = 0; i < dsc->var_cnt; i++)
 904     {
 905       struct dsc_var *dv = &dsc->vars[i];
 906       size_t j;
 907
 908       nc = 0;
 909       tab_text (t, nc++, i + 1, TAB_LEFT, dv->v->name);
 910       tab_text (t, nc++, i + 1, TAT_PRINTF, "%g", dv->valid);
 911       if (dsc->format == DSC_SERIAL)
 912         tab_text (t, nc++, i + 1, TAT_PRINTF, "%g", dv->missing);
 913       for (j = 0; j < DSC_N_STATS; j++)
 914         if (dsc->show_stats & (1ul << j))
 915           tab_float (t, nc++, i + 1, TAB_NONE, dv->stats[j], 10, 3);
 916     }
 917
 918   tab_title (t, 1, _("Valid cases = %g; cases with missing value(s) = %g."),
 919              dsc->valid, dsc->missing_listwise);
 920
 921   tab_submit (t);
 922 }
 923
 924 /* Compares `struct dsc_var's A and B according to the ordering
 925    specified by CMD. */
 926 static int
 927 descriptives_compare_dsc_vars (const void *a_, const void *b_, void *dsc_)
 928 {
 929   const struct dsc_var *a = a_;
 930   const struct dsc_var *b = b_;
 931   struct dsc_proc *dsc = dsc_;
 932
 933   int result;
 934
 935   if (dsc->sort_by_stat == DSC_NAME)
 936     result = strcasecmp (a->v->name, b->v->name);
 937   else
 938     {
 939       double as = a->stats[dsc->sort_by_stat];
 940       double bs = b->stats[dsc->sort_by_stat];
 941
 942       result = as < bs ? -1 : as > bs;
 943     }
 944
 945   if (!dsc->sort_ascending)
 946     result = -result;
 947
 948   return result;
 949 }