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