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