src/language/stats/descriptives.c

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