src/descript.q

   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 <assert.h>
  24 #include <limits.h>
  25 #include <math.h>
  26 #include <stdlib.h>
  27 #include "algorithm.h"
  28 #include "alloc.h"
  29 #include "bitvector.h"
  30 #include "command.h"
  31 #include "lexer.h"
  32 #include "error.h"
  33 #include "approx.h"
  34 #include "magic.h"
  35 #include "stats.h"
  36 #include "som.h"
  37 #include "tab.h"
  38 #include "var.h"
  39 #include "vfm.h"
  40
  41 /* (specification)
  42    DESCRIPTIVES (dsc_):
  43      *variables=custom;
  44      +missing=miss:!variable/listwise,incl:!noinclude/include;
  45      +format=labeled:!labels/nolabels,indexed:!noindex/index,lined:!line/serial;
  46      +save=;
  47      +options[op_]=1,2,3,4,5,6,7,8;
  48      +statistics[st_]=all,1|mean,2|semean,5|stddev,6|variance,7|kurtosis,
  49                       8|skewness,9|range,10|minimum,11|maximum,12|sum,
  50                       13|default,14|seskewness,15|sekurtosis;
  51      +sort=sortby:mean/semean/stddev/variance/kurtosis/skewness/range/
  52            range/minimum/maximum/sum/name/seskewness/sekurtosis/!none,
  53            order:!a/d.
  54 */
  55 /* (declarations) */
  56 /* (functions) */
  57
  58 /* DESCRIPTIVES private data. */
  59
  60 /* Describes properties of a distribution for the purpose of
  61    calculating a Z-score. */
  62 struct dsc_z_score
  63   {
  64     struct variable *s, *d;     /* Source, destination variable. */
  65     double mean;                /* Distribution mean. */
  66     double std_dev;             /* Distribution standard deviation. */
  67   };
  68
  69 /* DESCRIPTIVES transformation (for calculating Z-scores). */
  70 struct descriptives_trns
  71   {
  72     struct trns_header h;
  73     int n;                      /* Number of Z-scores. */
  74     struct dsc_z_score *z;      /* Array of Z-scores. */
  75   };
  76
  77 /* These next three vars, see comment at top of display(). */
  78 /* Number of cases missing listwise, even if option 5 not selected. */
  79 static double d_glob_miss_list;
  80
  81 /* Number of total *cases* valid or missing, as a double.  Unless
  82    option 5 is selected, d_glob_missing is 0. */
  83 static double d_glob_valid, d_glob_missing;
  84
  85 /* Set when a weighting variable is missing or <=0. */
  86 static int bad_weight;
  87
  88 /* Number of generic zvarnames we've generated in this execution. */
  89 static int z_count;
  90
  91 /* Variables specified on command. */
  92 static struct variable **v_variables;
  93 static int n_variables;
  94
  95 /* Command specifications. */
  96 static struct cmd_descriptives cmd;
  97
  98 /* Whether z-scores are computed. */
  99 static int z_scores;
 100
 101 /* Statistic to sort by. */
 102 static int sortby_stat;
 103
 104 /* Statistics to display. */
 105 static unsigned long stats;
 106
 107 /* Easier access to long-named arrays. */
 108 #define stat cmd.a_statistics
 109 #define opt  cmd.a_options
 110
 111 /* Groups of statistics. */
 112 #define BI          BIT_INDEX
 113
 114 #define dsc_default                                                     \
 115         (BI (dsc_mean) | BI (dsc_stddev) | BI (dsc_min) | BI (dsc_max))
 116
 117 #define dsc_all                                                 \
 118         (BI (dsc_sum) | BI (dsc_min) | BI (dsc_max)             \
 119          | BI (dsc_mean) | BI (dsc_semean) | BI (dsc_stddev)    \
 120          | BI (dsc_variance) | BI (dsc_kurt) | BI (dsc_sekurt)  \
 121          | BI (dsc_skew) | BI (dsc_seskew) | BI (dsc_range)     \
 122          | BI (dsc_range))
 123
 124 /* Table of options. */
 125 #define op_incl_miss    DSC_OP_1        /* Honored. */
 126 #define op_no_varlabs   DSC_OP_2        /* Ignored. */
 127 #define op_zscores      DSC_OP_3        /* Honored. */
 128 #define op_index        DSC_OP_4        /* FIXME. */
 129 #define op_excl_miss    DSC_OP_5        /* Honored. */
 130 #define op_serial       DSC_OP_6        /* Honored. */
 131 #define op_narrow       DSC_OP_7        /* Ignored. */
 132 #define op_no_varnames  DSC_OP_8        /* Honored. */
 133
 134 /* Describes one statistic that can be calculated. */
 135 /* FIXME: Currently sm,col_width are not used. */
 136 struct dsc_info
 137   {
 138     int st_indx;                /* Index into st_a_statistics[]. */
 139     int sb_indx;                /* Sort-by index. */
 140     const char *s10;            /* Name, stuffed into 10 columns. */
 141     const char *s8;             /* Name, stuffed into 8 columns. */
 142     const char *sm;             /* Name, stuffed minimally. */
 143     const char *s;              /* Full name. */
 144     int max_degree;             /* Highest degree necessary to calculate this
 145                                    statistic. */
 146     int col_width;              /* Column width (not incl. spacing between columns) */
 147   };
 148
 149 /* Table of statistics, indexed by dsc_*. */
 150 static struct dsc_info dsc_info[dsc_n_stats] =
 151 {
 152   {DSC_ST_MEAN, DSC_MEAN, N_("Mean"), N_("Mean"), N_("Mean"), N_("mean"), 1, 10},
 153   {DSC_ST_SEMEAN, DSC_SEMEAN, N_("S.E. Mean"), N_("S E Mean"), N_("SE"),
 154    N_("standard error of mean"), 2, 10},
 155   {DSC_ST_STDDEV, DSC_STDDEV, N_("Std Dev"), N_("Std Dev"), N_("SD"),
 156    N_("standard deviation"), 2, 11},
 157   {DSC_ST_VARIANCE, DSC_VARIANCE, N_("Variance"), N_("Variance"),
 158    N_("Var"), N_("variance"), 2, 12},
 159   {DSC_ST_KURTOSIS, DSC_KURTOSIS, N_("Kurtosis"), N_("Kurtosis"),
 160    N_("Kurt"), N_("kurtosis"), 4, 9},
 161   {DSC_ST_SEKURTOSIS, DSC_SEKURTOSIS, N_("S.E. Kurt"), N_("S E Kurt"), N_("SEKurt"),
 162    N_("standard error of kurtosis"), 0, 9},
 163   {DSC_ST_SKEWNESS, DSC_SKEWNESS, N_("Skewness"), N_("Skewness"), N_("Skew"),
 164    N_("skewness"), 3, 9},
 165   {DSC_ST_SESKEWNESS, DSC_SESKEWNESS, N_("S.E. Skew"), N_("S E Skew"), N_("SESkew"),
 166    N_("standard error of skewness"), 0, 9},
 167   {DSC_ST_RANGE, DSC_RANGE, N_("Range"), N_("Range"), N_("Rng"), N_("range"), 0, 10},
 168   {DSC_ST_MINIMUM, DSC_MINIMUM, N_("Minimum"), N_("Minimum"), N_("Min"),
 169    N_("minimum"), 0, 10},
 170   {DSC_ST_MAXIMUM, DSC_MAXIMUM, N_("Maximum"), N_("Maximum"), N_("Max"),
 171    N_("maximum"), 0, 10},
 172   {DSC_ST_SUM, DSC_SUM, N_("Sum"), N_("Sum"), N_("Sum"), N_("sum"), 1, 13},
 173 };
 174
 175 /* Z-score functions. */
 176 static int generate_z_varname (struct variable * v);
 177 static void dump_z_table (void);
 178 static void run_z_pass (void);
 179
 180 /* Procedure execution functions. */
 181 static int calc (struct ccase *);
 182 static void precalc (void);
 183 static void postcalc (void);
 184 static void display (void);
 185 \f
 186 /* Parser and outline. */
 187
 188 int
 189 cmd_descriptives (void)
 190 {
 191   struct variable *v;
 192   int i;
 193
 194   v_variables = NULL;
 195   n_variables = 0;
 196
 197   lex_match_id ("DESCRIPTIVES");
 198   lex_match_id ("CONDESCRIPTIVES");
 199   if (!parse_descriptives (&cmd))
 200     goto lossage;
 201
 202   if (n_variables == 0)
 203     goto lossage;
 204   for (i = 0; i < n_variables; i++)
 205     {
 206       v = v_variables[i];
 207       v->p.dsc.dup = 0;
 208       v->p.dsc.zname[0] = 0;
 209     }
 210
 211   if (n_variables < 0)
 212     {
 213       msg (SE, _("No variables specified."));
 214       goto lossage;
 215     }
 216
 217   if (cmd.sbc_options && (cmd.sbc_save || cmd.sbc_format || cmd.sbc_missing))
 218     {
 219       msg (SE, _("OPTIONS may not be used with SAVE, FORMAT, or MISSING."));
 220       goto lossage;
 221     }
 222
 223   if (!cmd.sbc_options)
 224     {
 225       if (cmd.incl == DSC_INCLUDE)
 226         opt[op_incl_miss] = 1;
 227       if (cmd.labeled == DSC_NOLABELS)
 228         opt[op_no_varlabs] = 1;
 229       if (cmd.sbc_save)
 230         opt[op_zscores] = 1;
 231       if (cmd.miss == DSC_LISTWISE)
 232         opt[op_excl_miss] = 1;
 233       if (cmd.lined == DSC_SERIAL)
 234         opt[op_serial] = 1;
 235     }
 236
 237   /* Construct z-score varnames, show translation table. */
 238   if (opt[op_zscores])
 239     {
 240       z_count = 0;
 241       for (i = 0; i < n_variables; i++)
 242         {
 243           v = v_variables[i];
 244           if (v->p.dsc.dup++)
 245             continue;
 246
 247           if (v->p.dsc.zname[0] == 0)
 248             if (!generate_z_varname (v))
 249               goto lossage;
 250         }
 251       dump_z_table ();
 252       z_scores = 1;
 253     }
 254
 255   /* Figure out statistics to calculate. */
 256   stats = 0;
 257   if (stat[DSC_ST_DEFAULT] || !cmd.sbc_statistics)
 258     stats |= dsc_default;
 259   if (stat[DSC_ST_ALL])
 260     stats |= dsc_all;
 261   for (i = 0; i < dsc_n_stats; i++)
 262     if (stat[dsc_info[i].st_indx])
 263       stats |= BIT_INDEX (i);
 264   if (stats & dsc_kurt)
 265     stats |= dsc_sekurt;
 266   if (stats & dsc_skew)
 267     stats |= dsc_seskew;
 268
 269   /* Check the sort order. */
 270   sortby_stat = -1;
 271   if (cmd.sortby == DSC_NONE)
 272     sortby_stat = -2;
 273   else if (cmd.sortby != DSC_NAME)
 274     {
 275       for (i = 0; i < n_variables; i++)
 276         if (dsc_info[i].sb_indx == cmd.sortby)
 277           {
 278             sortby_stat = i;
 279             if (!(stats & BIT_INDEX (i)))
 280               {
 281                 msg (SE, _("It's not possible to sort on `%s' without "
 282                            "displaying `%s'."),
 283                      gettext (dsc_info[i].s), gettext (dsc_info[i].s));
 284                 goto lossage;
 285               }
 286           }
 287       assert (sortby_stat != -1);
 288     }
 289
 290   /* Data pass! */
 291   update_weighting (&default_dict);
 292   bad_weight = 0;
 293   procedure (precalc, calc, postcalc);
 294
 295   if (bad_weight)
 296     msg (SW, _("At least one case in the data file had a weight value "
 297          "that was system-missing, zero, or negative.  These case(s) "
 298          "were ignored."));
 299
 300   /* Z-scoring! */
 301   if (z_scores)
 302     run_z_pass ();
 303
 304   if (v_variables)
 305     free (v_variables);
 306   return CMD_SUCCESS;
 307
 308  lossage:
 309   if (v_variables)
 310     free (v_variables);
 311   return CMD_FAILURE;
 312 }
 313
 314 /* Parses the VARIABLES subcommand. */
 315 static int
 316 dsc_custom_variables (struct cmd_descriptives *cmd unused)
 317 {
 318   if (!lex_match_id ("VARIABLES")
 319       && (token != T_ID || !is_varname (tokid))
 320       && token != T_ALL)
 321     return 2;
 322   lex_match ('=');
 323
 324   while (token == T_ID || token == T_ALL)
 325     {
 326       int i, n;
 327
 328       n = n_variables;
 329       if (!parse_variables (NULL, &v_variables, &n_variables,
 330                             PV_DUPLICATE | PV_SINGLE | PV_APPEND | PV_NUMERIC
 331                             | PV_NO_SCRATCH))
 332         return 0;
 333       if (lex_match ('('))
 334         {
 335           if (n_variables - n > 1)
 336             {
 337               msg (SE, _("Names for z-score variables must be given for "
 338                          "individual variables, not for groups of "
 339                          "variables."));
 340               return 0;
 341             }
 342           assert (n_variables - n <= 0);
 343           if (token != T_ID)
 344             {
 345               msg (SE, _("Name for z-score variable expected."));
 346               return 0;
 347             }
 348           if (is_varname (tokid))
 349             {
 350               msg (SE, _("Z-score variable name `%s' is a "
 351                          "duplicate variable name with a current variable."),
 352                    tokid);
 353               return 0;
 354             }
 355           for (i = 0; i < n_variables; i++)
 356             if (v_variables[i]->p.dsc.zname[0]
 357                 && !strcmp (v_variables[i]->p.dsc.zname, tokid))
 358               {
 359                 msg (SE, _("Z-score variable name `%s' is "
 360                            "used multiple times."), tokid);
 361                 return 0;
 362               }
 363           strcpy (v_variables[n_variables - 1]->p.dsc.zname, tokid);
 364           lex_get ();
 365           if (token != ')')
 366             {
 367               msg (SE, _("`)' expected after z-score variable name."));
 368               return 0;
 369             }
 370
 371           z_scores = 1;
 372         }
 373       lex_match (',');
 374     }
 375   return 1;
 376 }
 377 \f
 378 /* Z scores. */
 379
 380 /* Returns 0 if NAME is a duplicate of any existing variable name or
 381    of any previously-declared z-var name; otherwise returns 1. */
 382 static int
 383 try_name (char *name)
 384 {
 385   int i;
 386
 387   if (is_varname (name))
 388     return 0;
 389   for (i = 0; i < n_variables; i++)
 390     {
 391       struct variable *v = v_variables[i];
 392       if (!strcmp (v->p.dsc.zname, name))
 393         return 0;
 394     }
 395   return 1;
 396 }
 397
 398 static int
 399 generate_z_varname (struct variable * v)
 400 {
 401   char zname[10];
 402
 403   strcpy (&zname[1], v->name);
 404   zname[0] = 'Z';
 405   zname[8] = '\0';
 406   if (try_name (zname))
 407     {
 408       strcpy (v->p.dsc.zname, zname);
 409       return 1;
 410     }
 411
 412   for (;;)
 413     {
 414       /* Generate variable name. */
 415       z_count++;
 416
 417       if (z_count <= 99)
 418         sprintf (zname, "ZSC%03d", z_count);
 419       else if (z_count <= 108)
 420         sprintf (zname, "STDZ%02d", z_count - 99);
 421       else if (z_count <= 117)
 422         sprintf (zname, "ZZZZ%02d", z_count - 108);
 423       else if (z_count <= 126)
 424         sprintf (zname, "ZQZQ%02d", z_count - 117);
 425       else
 426         {
 427           msg (SE, _("Ran out of generic names for Z-score variables.  "
 428                      "There are only 126 generic names: ZSC001-ZSC0999, "
 429                      "STDZ01-STDZ09, ZZZZ01-ZZZZ09, ZQZQ01-ZQZQ09."));
 430           return 0;
 431         }
 432
 433       if (try_name (zname))
 434         {
 435           strcpy (v->p.dsc.zname, zname);
 436           return 1;
 437         }
 438     }
 439 }
 440
 441 static void
 442 dump_z_table (void)
 443 {
 444   int count;
 445   struct tab_table *t;
 446
 447   {
 448     int i;
 449
 450     for (count = i = 0; i < n_variables; i++)
 451       if (v_variables[i]->p.dsc.zname)
 452         count++;
 453   }
 454
 455   t = tab_create (2, count + 1, 0);
 456   tab_title (t, 0, _("Mapping of variables to corresponding Z-scores."));
 457   tab_columns (t, SOM_COL_DOWN, 1);
 458   tab_headers (t, 0, 0, 1, 0);
 459   tab_box (t, TAL_1, TAL_1, TAL_0, TAL_1, 0, 0, 1, count);
 460   tab_hline (t, TAL_2, 0, 1, 1);
 461   tab_text (t, 0, 0, TAB_CENTER | TAT_TITLE, _("Source"));
 462   tab_text (t, 1, 0, TAB_CENTER | TAT_TITLE, _("Target"));
 463   tab_dim (t, tab_natural_dimensions);
 464
 465   {
 466     int i, y;
 467
 468     for (i = 0, y = 1; i < n_variables; i++)
 469       if (v_variables[i]->p.dsc.zname)
 470         {
 471           tab_text (t, 0, y, TAB_LEFT, v_variables[i]->name);
 472           tab_text (t, 1, y++, TAB_LEFT, v_variables[i]->p.dsc.zname);
 473         }
 474   }
 475
 476   tab_submit (t);
 477 }
 478
 479 /* Transformation function to calculate Z-scores. */
 480 static int
 481 descriptives_trns_proc (struct trns_header * trns, struct ccase * c)
 482 {
 483   struct descriptives_trns *t = (struct descriptives_trns *) trns;
 484   struct dsc_z_score *z;
 485   int i;
 486
 487   for (i = t->n, z = t->z; i--; z++)
 488     {
 489       double score = c->data[z->s->fv].f;
 490
 491       if (z->mean == SYSMIS || score == SYSMIS)
 492         c->data[z->d->fv].f = SYSMIS;
 493       else
 494         c->data[z->d->fv].f = (score - z->mean) / z->std_dev;
 495     }
 496   return -1;
 497 }
 498
 499 /* Frees a descriptives_trns struct. */
 500 static void
 501 descriptives_trns_free (struct trns_header * trns)
 502 {
 503   struct descriptives_trns *t = (struct descriptives_trns *) trns;
 504
 505   free (t->z);
 506 }
 507
 508 /* The name is a misnomer: actually this function sets up a
 509    transformation by which scores can be converted into Z-scores. */
 510 static void
 511 run_z_pass (void)
 512 {
 513   struct descriptives_trns *t;
 514   int count, i;
 515
 516   for (i = 0; i < n_variables; i++)
 517     v_variables[i]->p.dsc.dup = 0;
 518   for (count = i = 0; i < n_variables; i++)
 519     {
 520       if (v_variables[i]->p.dsc.dup++)
 521         continue;
 522       if (v_variables[i]->p.dsc.zname)
 523         count++;
 524     }
 525
 526   t = xmalloc (sizeof *t);
 527   t->h.proc = descriptives_trns_proc;
 528   t->h.free = descriptives_trns_free;
 529   t->n = count;
 530   t->z = xmalloc (count * sizeof *t->z);
 531
 532   for (i = 0; i < n_variables; i++)
 533     v_variables[i]->p.dsc.dup = 0;
 534   for (count = i = 0; i < n_variables; i++)
 535     {
 536       struct variable *v = v_variables[i];
 537       if (v->p.dsc.dup++ == 0 && v->p.dsc.zname[0])
 538         {
 539           char *cp;
 540           struct variable *d;
 541
 542           t->z[count].s = v;
 543           t->z[count].d = d = force_create_variable (&default_dict,
 544                                                      v->p.dsc.zname,
 545                                                      NUMERIC, 0);
 546           if (v->label)
 547             {
 548               d->label = xmalloc (strlen (v->label) + 12);
 549               cp = stpcpy (d->label, _("Z-score of "));
 550               strcpy (cp, v->label);
 551             }
 552           else
 553             {
 554               d->label = xmalloc (strlen (v->name) + 12);
 555               cp = stpcpy (d->label, _("Z-score of "));
 556               strcpy (cp, v->name);
 557             }
 558           t->z[count].mean = v->p.dsc.stats[dsc_mean];
 559           t->z[count].std_dev = v->p.dsc.stats[dsc_stddev];
 560           if (t->z[count].std_dev == SYSMIS
 561               || approx_eq (t->z[count].std_dev, 0.0))
 562             t->z[count].mean = SYSMIS;
 563           count++;
 564         }
 565     }
 566
 567   add_transformation ((struct trns_header *) t);
 568 }
 569 \f
 570 /* Statistical calculation. */
 571
 572 static void
 573 precalc (void)
 574 {
 575   int i;
 576
 577   for (i = 0; i < n_variables; i++)
 578     v_variables[i]->p.dsc.dup = -2;
 579   for (i = 0; i < n_variables; i++)
 580     {
 581       struct descriptives_proc *dsc = &v_variables[i]->p.dsc;
 582
 583       /* Don't need to initialize more than once. */
 584       if (dsc->dup == -1)
 585         continue;
 586       dsc->dup = -1;
 587
 588       dsc->valid = dsc->miss = 0.0;
 589       dsc->X_bar = dsc->M2 = dsc->M3 = dsc->M4 = 0.0;
 590       dsc->min = DBL_MAX;
 591       dsc->max = -DBL_MAX;
 592     }
 593   d_glob_valid = d_glob_missing = 0.0;
 594 }
 595
 596 static int
 597 calc (struct ccase * c)
 598 {
 599   int i;
 600
 601   /* Unique case identifier. */
 602   static int case_id;
 603
 604   /* Get the weight for this case. */
 605   double w;
 606
 607   if (default_dict.weight_index == -1)
 608     w = 1.0;
 609   else
 610     {
 611       w = c->data[default_dict.weight_index].f;
 612       if (w <= 0.0 || w == SYSMIS)
 613         {
 614           w = 0.0;
 615           bad_weight = 1;
 616         }
 617     }
 618
 619   case_id++;
 620
 621   /* Handle missing values. */
 622   for (i = 0; i < n_variables; i++)
 623     {
 624       struct variable *v = v_variables[i];
 625       double X = c->data[v->fv].f;
 626
 627       if (X == SYSMIS || (!opt[op_incl_miss] && is_num_user_missing (X, v)))
 628         {
 629           if (opt[op_excl_miss])
 630             {
 631               d_glob_missing += w;
 632               return 1;
 633             }
 634           else
 635             {
 636               d_glob_miss_list += w;
 637               goto iterate;
 638             }
 639         }
 640     }
 641   d_glob_valid += w;
 642
 643 iterate:
 644   for (i = 0; i < n_variables; i++)
 645     {
 646       struct descriptives_proc *inf = &v_variables[i]->p.dsc;
 647
 648       double X, v;
 649       double W_old, W_new;
 650       double v2, v3, v4;
 651       double w2, w3, w4;
 652
 653       if (inf->dup == case_id)
 654         continue;
 655       inf->dup = case_id;
 656
 657       X = c->data[v_variables[i]->fv].f;
 658       if (X == SYSMIS
 659           || (!opt[op_incl_miss] && is_num_user_missing (X, v_variables[i])))
 660         {
 661           inf->miss += w;
 662           continue;
 663         }
 664
 665       /* These formulas taken from _SPSS Statistical Algorithms_.  The
 666          names W_old, and W_new are used for W_j-1 and W_j,
 667          respectively, and other variables simply have the subscripts
 668          trimmed off, except for X_bar.
 669
 670          I am happy that mathematical formulas may not be
 671          copyrighted. */
 672       W_old = inf->valid;
 673       W_new = inf->valid += w;
 674       v = (w / W_new) * (X - inf->X_bar);
 675       v2 = v * v;
 676       v3 = v2 * v;
 677       v4 = v3 * v;
 678       w2 = w * w;
 679       w3 = w2 * w;
 680       w4 = w3 * w;
 681       inf->M4 += (-4.0 * v * inf->M3 + 6.0 * v2 * inf->M2
 682                + (W_new * W_new - 3 * w * W_old / w3) * v4 * W_old * W_new);
 683       inf->M3 += (-3.0 * v * inf->M2 + W_new * W_old / w2
 684                   * (W_new - 2 * w) * v3);
 685       inf->M2 += W_new * W_old / w * v2;
 686       inf->X_bar += v;
 687       if (X < inf->min)
 688         inf->min = X;
 689       if (X > inf->max)
 690         inf->max = X;
 691     }
 692   return 1;
 693 }
 694
 695 static void
 696 postcalc (void)
 697 {
 698   int i;
 699
 700   if (opt[op_excl_miss])
 701     d_glob_miss_list = d_glob_missing;
 702
 703   for (i = 0; i < n_variables; i++)
 704     {
 705       struct descriptives_proc *dsc = &v_variables[i]->p.dsc;
 706       double W;
 707
 708       /* Don't duplicate our efforts. */
 709       if (dsc->dup == -2)
 710         continue;
 711       dsc->dup = -2;
 712
 713       W = dsc->valid;
 714
 715       dsc->stats[dsc_mean] = dsc->X_bar;
 716       dsc->stats[dsc_variance] = dsc->M2 / (W - 1);
 717       dsc->stats[dsc_stddev] = sqrt (dsc->stats[dsc_variance]);
 718       dsc->stats[dsc_semean] = dsc->stats[dsc_stddev] / sqrt (W);
 719       dsc->stats[dsc_min] = dsc->min == DBL_MAX ? SYSMIS : dsc->min;
 720       dsc->stats[dsc_max] = dsc->max == -DBL_MAX ? SYSMIS : dsc->max;
 721       dsc->stats[dsc_range] = ((dsc->min == DBL_MAX || dsc->max == -DBL_MAX)
 722                                ? SYSMIS : dsc->max - dsc->min);
 723       dsc->stats[dsc_sum] = W * dsc->X_bar;
 724       if (W > 2.0 && dsc->stats[dsc_variance] >= 1e-20)
 725         {
 726           double S = dsc->stats[dsc_stddev];
 727           dsc->stats[dsc_skew] = (W * dsc->M3 / ((W - 1.0) * (W - 2.0) * S * S * S));
 728           dsc->stats[dsc_seskew] =
 729             sqrt (6.0 * W * (W - 1.0) / ((W - 2.0) * (W + 1.0) * (W + 3.0)));
 730         }
 731       else
 732         {
 733           dsc->stats[dsc_skew] = dsc->stats[dsc_seskew] = SYSMIS;
 734         }
 735       if (W > 3.0 && dsc->stats[dsc_variance] >= 1e-20)
 736         {
 737           double S2 = dsc->stats[dsc_variance];
 738           double SE_g1 = dsc->stats[dsc_seskew];
 739
 740           dsc->stats[dsc_kurt] =
 741             (W * (W + 1.0) * dsc->M4 - 3.0 * dsc->M2 * dsc->M2 * (W - 1.0))
 742             / ((W - 1.0) * (W - 2.0) * (W - 3.0) * S2 * S2);
 743
 744           /* Note that in _SPSS Statistical Algorithms_, the square
 745              root symbol is missing from this formula. */
 746           dsc->stats[dsc_sekurt] =
 747             sqrt ((4.0 * (W * W - 1.0) * SE_g1 * SE_g1) / ((W - 3.0) * (W + 5.0)));
 748         }
 749       else
 750         {
 751           dsc->stats[dsc_kurt] = dsc->stats[dsc_sekurt] = SYSMIS;
 752         }
 753     }
 754
 755   display ();
 756 }
 757 \f
 758 /* Statistical display. */
 759
 760 static algo_compare_func descriptives_compare_variables;
 761
 762 static void
 763 display (void)
 764 {
 765   int i, j;
 766
 767   int nc, n_stats;
 768   struct tab_table *t;
 769
 770   /* If op_excl_miss is on, d_glob_valid and (potentially)
 771      d_glob_missing are nonzero, and d_glob_missing equals
 772      d_glob_miss_list.
 773
 774      If op_excl_miss is off, d_glob_valid is nonzero.  d_glob_missing
 775      is zero.  d_glob_miss_list is (potentially) nonzero.  */
 776
 777   if (sortby_stat != -2)
 778     sort (v_variables, n_variables, sizeof *v_variables,
 779            descriptives_compare_variables, &cmd);
 780
 781   for (nc = i = 0; i < dsc_n_stats; i++)
 782     if (stats & BIT_INDEX (i))
 783       nc++;
 784   n_stats = nc;
 785   if (!opt[op_no_varnames])
 786     nc++;
 787   nc += opt[op_serial] ? 2 : 1;
 788
 789   t = tab_create (nc, n_variables + 1, 0);
 790   tab_headers (t, 1, 0, 1, 0);
 791   tab_box (t, TAL_1, TAL_1, -1, -1, 0, 0, nc - 1, n_variables);
 792   tab_box (t, -1, -1, -1, TAL_1, 1, 0, nc - 1, n_variables);
 793   tab_hline (t, TAL_2, 0, nc - 1, 1);
 794   tab_vline (t, TAL_2, 1, 0, n_variables);
 795   tab_dim (t, tab_natural_dimensions);
 796
 797   nc = 0;
 798   if (!opt[op_no_varnames])
 799     {
 800       tab_text (t, nc++, 0, TAB_LEFT | TAT_TITLE, _("Variable"));
 801     }
 802   if (opt[op_serial])
 803     {
 804       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Valid N"));
 805       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, _("Missing N"));
 806     } else {
 807       tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, "N");
 808     }
 809
 810   for (i = 0; i < dsc_n_stats; i++)
 811     if (stats & BIT_INDEX (i))
 812       {
 813         const char *title = gettext (dsc_info[i].s8);
 814         tab_text (t, nc++, 0, TAB_CENTER | TAT_TITLE, title);
 815       }
 816
 817   for (i = 0; i < n_variables; i++)
 818     {
 819       struct variable *v = v_variables[i];
 820
 821       nc = 0;
 822       if (!opt[op_no_varnames])
 823         tab_text (t, nc++, i + 1, TAB_LEFT, v->name);
 824       tab_text (t, nc++, i + 1, TAT_PRINTF, "%g", v->p.dsc.valid);
 825       if (opt[op_serial])
 826         tab_text (t, nc++, i + 1, TAT_PRINTF, "%g", v->p.dsc.miss);
 827       for (j = 0; j < dsc_n_stats; j++)
 828         if (stats & BIT_INDEX (j))
 829           tab_float (t, nc++, i + 1, TAB_NONE, v->p.dsc.stats[j], 10, 3);
 830     }
 831
 832   tab_title (t, 1, _("Valid cases = %g; cases with missing value(s) = %g."),
 833              d_glob_valid, d_glob_miss_list);
 834
 835   tab_submit (t);
 836 }
 837
 838 /* Compares variables A and B according to the ordering specified
 839    by CMD. */
 840 static int
 841 descriptives_compare_variables (const void *a_, const void *b_, void *cmd_)
 842 {
 843   struct variable *const *ap = a_;
 844   struct variable *const *bp = b_;
 845   struct variable *a = *ap;
 846   struct variable *b = *bp;
 847   struct cmd_descriptives *cmd = cmd_;
 848
 849   int result;
 850
 851   if (cmd->sortby == DSC_NAME)
 852     result = strcmp (a->name, b->name);
 853   else
 854     {
 855       double as = a->p.dsc.stats[sortby_stat];
 856       double bs = b->p.dsc.stats[sortby_stat];
 857
 858       result = as < bs ? -1 : as > bs;
 859     }
 860
 861   if (cmd->order == DSC_D)
 862     result = -result;
 863
 864   return result;
 865 }
 866
 867 /*
 868    Local variables:
 869    mode: c
 870    End:
 871 */