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