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