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