src/moments.c

   1 /* PSPP - computes sample statistics.
   2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
   3    Written by Ben Pfaff <blp@gnu.org>.
   4
   5    This program is free software; you can redistribute it and/or
   6    modify it under the terms of the GNU General Public License as
   7    published by the Free Software Foundation; either version 2 of the
   8    License, or (at your option) any later version.
   9
  10    This program is distributed in the hope that it will be useful, but
  11    WITHOUT ANY WARRANTY; without even the implied warranty of
  12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  13    General Public License for more details.
  14
  15    You should have received a copy of the GNU General Public License
  16    along with this program; if not, write to the Free Software
  17    Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
  18    02110-1301, USA. */
  19
  20 #include <config.h>
  21 #include "moments.h"
  22 #include <assert.h>
  23 #include <math.h>
  24 #include <stdlib.h>
  25 #include "alloc.h"
  26 #include "misc.h"
  27 #include "val.h"
  28
  29 #include "gettext.h"
  30 #define _(msgid) gettext (msgid)
  31 \f
  32 /* Calculates variance, skewness, and kurtosis into *VARIANCE,
  33    *SKEWNESS, and *KURTOSIS if they are non-null and not greater
  34    moments than MAX_MOMENT.  Accepts W as the total weight, D1 as
  35    the total deviation from the estimated mean, and D2, D3, and
  36    D4 as the sum of the squares, cubes, and 4th powers,
  37    respectively, of the deviation from the estimated mean. */
  38 static void
  39 calc_moments (enum moment max_moment,
  40               double w, double d1, double d2, double d3, double d4,
  41               double *variance, double *skewness, double *kurtosis)
  42 {
  43   assert (w > 0.);
  44
  45   if (max_moment >= MOMENT_VARIANCE && w > 1.)
  46     {
  47       double s2;
  48
  49       /* From _Numerical Recipes in C_, 2nd ed., 0-521-43108-5,
  50          section 14.1. */
  51       s2 = (d2 - pow2 (d1) / w) / (w - 1.);
  52       if (variance != NULL)
  53         *variance = s2;
  54
  55       /* From _SPSS Statistical Algorithms, 2nd ed.,
  56          0-918469-89-9, section "DESCRIPTIVES". */
  57       if (fabs (*variance) >= 1e-20)
  58         {
  59           if (max_moment >= MOMENT_SKEWNESS && skewness != NULL && w > 2.)
  60             {
  61               double s3 = s2 * sqrt (s2);
  62               double g1 = (w * d3) / ((w - 1.0) * (w - 2.0) * s3);
  63               if (finite (g1))
  64                 *skewness = g1;
  65             }
  66           if (max_moment >= MOMENT_KURTOSIS && kurtosis != NULL && w > 3.)
  67             {
  68               double den = (w - 2.) * (w - 3.) * pow2 (s2);
  69               double g2 = (w * (w + 1) * d4 / (w - 1.) / den
  70                            - 3. * pow2 (d2) / den);
  71               if (finite (g2))
  72                 *kurtosis = g2;
  73             }
  74         }
  75     }
  76 }
  77 \f
  78 /* Two-pass moments. */
  79
  80 /* A set of two-pass moments. */
  81 struct moments
  82   {
  83     enum moment max_moment;     /* Highest-order moment we're computing. */
  84     int pass;                   /* Current pass (1 or 2). */
  85
  86     /* Pass one. */
  87     double w1;                  /* Total weight for pass 1, so far. */
  88     double sum;                 /* Sum of values so far. */
  89     double mean;                /* Mean = sum / w1. */
  90
  91     /* Pass two. */
  92     double w2;                  /* Total weight for pass 2, so far. */
  93     double d1;                  /* Sum of deviations from the mean. */
  94     double d2;                  /* Sum of squared deviations from the mean. */
  95     double d3;                  /* Sum of cubed deviations from the mean. */
  96     double d4;                  /* Sum of (deviations from the mean)**4. */
  97   };
  98
  99 /* Initializes moments M for calculating moment MAX_MOMENT and
 100    lower moments. */
 101 static void
 102 init_moments (struct moments *m, enum moment max_moment)
 103 {
 104   assert (m != NULL);
 105   assert (max_moment == MOMENT_MEAN || max_moment == MOMENT_VARIANCE
 106           || max_moment == MOMENT_SKEWNESS || max_moment == MOMENT_KURTOSIS);
 107   m->max_moment = max_moment;
 108   moments_clear (m);
 109 }
 110
 111 /* Clears out a set of moments so that it can be reused for a new
 112    set of values.  The moments to be calculated are not changed. */
 113 void
 114 moments_clear (struct moments *m)
 115 {
 116   m->pass = 1;
 117   m->w1 = m->w2 = 0.;
 118   m->sum = 0.;
 119 }
 120
 121 /* Creates and returns a data structure for calculating moment
 122    MAX_MOMENT and lower moments on a data series.  The user
 123    should call moments_pass_one() for each value in the series,
 124    then call moments_pass_two() for the same set of values in the
 125    same order, then call moments_calculate() to obtain the
 126    moments.  The user may ask for the mean at any time during the
 127    first pass (using moments_calculate()), but otherwise no
 128    statistics may be requested until the end of the second pass.
 129    Call moments_destroy() when the moments are no longer
 130    needed. */
 131 struct moments *
 132 moments_create (enum moment max_moment)
 133 {
 134   struct moments *m = xmalloc (sizeof *m);
 135   init_moments (m, max_moment);
 136   return m;
 137 }
 138
 139 /* Adds VALUE with the given WEIGHT to the calculation of
 140    moments for the first pass. */
 141 void
 142 moments_pass_one (struct moments *m, double value, double weight)
 143 {
 144   assert (m != NULL);
 145   assert (m->pass == 1);
 146
 147   if (value != SYSMIS && weight > 0.)
 148     {
 149       m->sum += value * weight;
 150       m->w1 += weight;
 151     }
 152 }
 153
 154 /* Adds VALUE with the given WEIGHT to the calculation of
 155    moments for the second pass. */
 156 void
 157 moments_pass_two (struct moments *m, double value, double weight)
 158 {
 159   double d, d_power;
 160
 161   assert (m != NULL);
 162
 163   if (m->pass == 1)
 164     {
 165       m->pass = 2;
 166       m->mean = m->w1 != 0. ? m->sum / m->w1 : 0.;
 167       m->d1 = m->d2 = m->d3 = m->d4 = 0.;
 168     }
 169
 170   if (value != SYSMIS && weight >= 0.)
 171     {
 172       m->w2 += weight;
 173
 174       d = d_power = value - m->mean;
 175       m->d1 += d_power * weight;
 176
 177       if (m->max_moment >= MOMENT_VARIANCE)
 178         {
 179           d_power *= d;
 180           m->d2 += d_power * weight;
 181
 182           if (m->max_moment >= MOMENT_SKEWNESS)
 183             {
 184               d_power *= d;
 185               m->d3 += d_power * weight;
 186
 187               if (m->max_moment >= MOMENT_KURTOSIS)
 188                 {
 189                   d_power *= d;
 190                   m->d4 += d_power * weight;
 191                 }
 192             }
 193         }
 194     }
 195 }
 196
 197 /* Calculates moments based on the input data.  Stores the total
 198    weight in *WEIGHT, the mean in *MEAN, the variance in
 199    *VARIANCE, the skewness in *SKEWNESS, and the kurtosis in
 200    *KURTOSIS.  Any of these result parameters may be null
 201    pointers, in which case the values are not calculated.  If any
 202    result cannot be calculated, either because they are undefined
 203    based on the input data or because their moments are higher
 204    than the maximum requested on moments_create(), then SYSMIS is
 205    stored into that result. */
 206 void
 207 moments_calculate (const struct moments *m,
 208                    double *weight,
 209                    double *mean, double *variance,
 210                    double *skewness, double *kurtosis)
 211 {
 212   assert (m != NULL);
 213
 214   if (mean != NULL)
 215     *mean = SYSMIS;
 216   if (variance != NULL)
 217     *variance = SYSMIS;
 218   if (skewness != NULL)
 219     *skewness = SYSMIS;
 220   if (kurtosis != NULL)
 221     *kurtosis = SYSMIS;
 222
 223   if (weight != NULL)
 224     *weight = m->w1;
 225
 226   /* How many passes so far? */
 227   if (m->pass == 1)
 228     {
 229       /* In the first pass we can only calculate the mean. */
 230       if (mean != NULL && m->w1 > 0.)
 231         *mean = m->sum / m->w1;
 232     }
 233   else
 234     {
 235       /* After the second pass we can calculate any stat.  We
 236          don't support "online" computation during the second
 237          pass, so As a simple self-check, the total weight for
 238          the passes must agree. */
 239       assert (m->pass == 2);
 240       assert (m->w1 == m->w2);
 241
 242       if (m->w2 > 0.)
 243         {
 244           if (mean != NULL)
 245             *mean = m->mean;
 246           calc_moments (m->max_moment,
 247                         m->w2, m->d1, m->d2, m->d3, m->d4,
 248                         variance, skewness, kurtosis);
 249         }
 250     }
 251 }
 252
 253 /* Destroys a set of moments. */
 254 void
 255 moments_destroy (struct moments *m)
 256 {
 257   free (m);
 258 }
 259
 260 /* Calculates the requested moments on the CNT values in ARRAY.
 261    Each value is given a weight of 1.  The total weight is stored
 262    into *WEIGHT (trivially) and the mean, variance, skewness, and
 263    kurtosis are stored into *MEAN, *VARIANCE, *SKEWNESS, and
 264    *KURTOSIS, respectively.  Any of the result pointers may be
 265    null, in which case no value is stored. */
 266 void
 267 moments_of_doubles (const double *array, size_t cnt,
 268                     double *weight,
 269                     double *mean, double *variance,
 270                     double *skewness, double *kurtosis)
 271 {
 272   enum moment max_moment;
 273   struct moments m;
 274   size_t idx;
 275
 276   if (kurtosis != NULL)
 277     max_moment = MOMENT_KURTOSIS;
 278   else if (skewness != NULL)
 279     max_moment = MOMENT_SKEWNESS;
 280   else if (variance != NULL)
 281     max_moment = MOMENT_VARIANCE;
 282   else
 283     max_moment = MOMENT_MEAN;
 284
 285   init_moments (&m, max_moment);
 286   for (idx = 0; idx < cnt; idx++)
 287     moments_pass_one (&m, array[idx], 1.);
 288   for (idx = 0; idx < cnt; idx++)
 289     moments_pass_two (&m, array[idx], 1.);
 290   moments_calculate (&m, weight, mean, variance, skewness, kurtosis);
 291 }
 292
 293 /* Calculates the requested moments on the CNT numeric values in
 294    ARRAY.  Each value is given a weight of 1.  The total weight
 295    is stored into *WEIGHT (trivially) and the mean, variance,
 296    skewness, and kurtosis are stored into *MEAN, *VARIANCE,
 297    *SKEWNESS, and *KURTOSIS, respectively.  Any of the result
 298    pointers may be null, in which case no value is stored. */
 299 void
 300 moments_of_values (const union value *array, size_t cnt,
 301                    double *weight,
 302                    double *mean, double *variance,
 303                    double *skewness, double *kurtosis)
 304 {
 305   enum moment max_moment;
 306   struct moments m;
 307   size_t idx;
 308
 309   if (kurtosis != NULL)
 310     max_moment = MOMENT_KURTOSIS;
 311   else if (skewness != NULL)
 312     max_moment = MOMENT_SKEWNESS;
 313   else if (variance != NULL)
 314     max_moment = MOMENT_VARIANCE;
 315   else
 316     max_moment = MOMENT_MEAN;
 317
 318   init_moments (&m, max_moment);
 319   for (idx = 0; idx < cnt; idx++)
 320     moments_pass_one (&m, array[idx].f, 1.);
 321   for (idx = 0; idx < cnt; idx++)
 322     moments_pass_two (&m, array[idx].f, 1.);
 323   moments_calculate (&m, weight, mean, variance, skewness, kurtosis);
 324 }
 325 \f
 326 /* One-pass moments. */
 327
 328 /* A set of one-pass moments. */
 329 struct moments1
 330   {
 331     enum moment max_moment;     /* Highest-order moment we're computing. */
 332     double w;                   /* Total weight so far. */
 333     double d1;                  /* Sum of deviations from the mean. */
 334     double d2;                  /* Sum of squared deviations from the mean. */
 335     double d3;                  /* Sum of cubed deviations from the mean. */
 336     double d4;                  /* Sum of (deviations from the mean)**4. */
 337   };
 338
 339 /* Initializes one-pass moments M for calculating moment
 340    MAX_MOMENT and lower moments. */
 341 static void
 342 init_moments1 (struct moments1 *m, enum moment max_moment)
 343 {
 344   assert (m != NULL);
 345   assert (max_moment == MOMENT_MEAN || max_moment == MOMENT_VARIANCE
 346           || max_moment == MOMENT_SKEWNESS || max_moment == MOMENT_KURTOSIS);
 347   m->max_moment = max_moment;
 348   moments1_clear (m);
 349 }
 350
 351 /* Clears out a set of one-pass moments so that it can be reused
 352    for a new set of values.  The moments to be calculated are not
 353    changed. */
 354 void
 355 moments1_clear (struct moments1 *m)
 356 {
 357   m->w = 0.;
 358   m->d1 = m->d2 = m->d3 = m->d4 = 0.;
 359 }
 360
 361 /* Creates and returns a data structure for calculating moment
 362    MAX_MOMENT and lower moments on a data series in a single
 363    pass.  The user should call moments1_add() for each value in
 364    the series.  The user may call moments1_calculate() to obtain
 365    the current moments at any time.  Call moments1_destroy() when
 366    the moments are no longer needed.
 367
 368    One-pass moments should only be used when two passes over the
 369    data are impractical. */
 370 struct moments1 *
 371 moments1_create (enum moment max_moment)
 372 {
 373   struct moments1 *m = xmalloc (sizeof *m);
 374   init_moments1 (m, max_moment);
 375   return m;
 376 }
 377
 378 /* Adds VALUE with the given WEIGHT to the calculation of
 379    one-pass moments. */
 380 void
 381 moments1_add (struct moments1 *m, double value, double weight)
 382 {
 383   assert (m != NULL);
 384
 385   if (value != SYSMIS && weight > 0.)
 386     {
 387       double prev_w, v1;
 388
 389       prev_w = m->w;
 390       m->w += weight;
 391       v1 = (weight / m->w) * (value - m->d1);
 392       m->d1 += v1;
 393
 394       if (m->max_moment >= MOMENT_VARIANCE)
 395         {
 396           double v2 = v1 * v1;
 397           double w_prev_w = m->w * prev_w;
 398           double prev_m2 = m->d2;
 399
 400           m->d2 += w_prev_w / weight * v2;
 401           if (m->max_moment >= MOMENT_SKEWNESS)
 402             {
 403               double w2 = weight * weight;
 404               double v3 = v2 * v1;
 405               double prev_m3 = m->d3;
 406
 407               m->d3 += (-3. * v1 * prev_m2
 408                          + w_prev_w / w2 * (m->w - 2. * weight) * v3);
 409               if (m->max_moment >= MOMENT_KURTOSIS)
 410                 {
 411                   double w3 = w2 * weight;
 412                   double v4 = v2 * v2;
 413
 414                   m->d4 += (-4. * v1 * prev_m3
 415                              + 6. * v2 * prev_m2
 416                              + ((pow2 (m->w) - 3. * weight * prev_w)
 417                                 * v4 * w_prev_w / w3));
 418                 }
 419             }
 420         }
 421     }
 422 }
 423
 424 /* Calculates one-pass moments based on the input data.  Stores
 425    the total weight in *WEIGHT, the mean in *MEAN, the variance
 426    in *VARIANCE, the skewness in *SKEWNESS, and the kurtosis in
 427    *KURTOSIS.  Any of these result parameters may be null
 428    pointers, in which case the values are not calculated.  If any
 429    result cannot be calculated, either because they are undefined
 430    based on the input data or because their moments are higher
 431    than the maximum requested on moments_create(), then SYSMIS is
 432    stored into that result. */
 433 void
 434 moments1_calculate (const struct moments1 *m,
 435                     double *weight,
 436                     double *mean, double *variance,
 437                     double *skewness, double *kurtosis)
 438 {
 439   assert (m != NULL);
 440
 441   if (mean != NULL)
 442     *mean = SYSMIS;
 443   if (variance != NULL)
 444     *variance = SYSMIS;
 445   if (skewness != NULL)
 446     *skewness = SYSMIS;
 447   if (kurtosis != NULL)
 448     *kurtosis = SYSMIS;
 449
 450   if (weight != NULL)
 451     *weight = m->w;
 452
 453   if (m->w > 0.)
 454     {
 455       if (mean != NULL)
 456         *mean = m->d1;
 457
 458       calc_moments (m->max_moment,
 459                     m->w, 0., m->d2, m->d3, m->d4,
 460                     variance, skewness, kurtosis);
 461     }
 462 }
 463
 464 /* Destroy one-pass moments M. */
 465 void
 466 moments1_destroy (struct moments1 *m)
 467 {
 468   free (m);
 469 }
 470 \f
 471 /* Returns the standard error of the skewness for the given total
 472    weight W.
 473
 474    From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
 475    section "DESCRIPTIVES". */
 476 double
 477 calc_seskew (double W)
 478 {
 479   return sqrt ((6. * W * (W - 1.)) / ((W - 2.) * (W + 1.) * (W + 3.)));
 480 }
 481
 482 /* Returns the standard error of the kurtosis for the given total
 483    weight W.
 484
 485    From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
 486    section "DESCRIPTIVES", except that the sqrt symbol is omitted
 487    there. */
 488 double
 489 calc_sekurt (double W)
 490 {
 491   return sqrt ((4. * (pow2 (W) - 1.) * pow2 (calc_seskew (W)))
 492                / ((W - 3.) * (W + 5.)));
 493 }
 494 \f
 495 #include <stdio.h>
 496 #include "command.h"
 497 #include "lexer.h"
 498
 499 static int
 500 read_values (double **values, double **weights, size_t *cnt)
 501 {
 502   size_t cap = 0;
 503
 504   *values = NULL;
 505   *weights = NULL;
 506   *cnt = 0;
 507   while (lex_is_number ())
 508     {
 509       double value = tokval;
 510       double weight = 1.;
 511       lex_get ();
 512       if (lex_match ('*'))
 513         {
 514           if (!lex_is_number ())
 515             {
 516               lex_error (_("expecting weight value"));
 517               return 0;
 518             }
 519           weight = tokval;
 520           lex_get ();
 521         }
 522
 523       if (*cnt >= cap)
 524         {
 525           cap = 2 * (cap + 8);
 526           *values = xnrealloc (*values, cap, sizeof **values);
 527           *weights = xnrealloc (*weights, cap, sizeof **weights);
 528         }
 529
 530       (*values)[*cnt] = value;
 531       (*weights)[*cnt] = weight;
 532       (*cnt)++;
 533     }
 534
 535   return 1;
 536 }
 537
 538 int
 539 cmd_debug_moments (void)
 540 {
 541   int retval = CMD_FAILURE;
 542   double *values = NULL;
 543   double *weights = NULL;
 544   double weight, M[4];
 545   int two_pass = 1;
 546   size_t cnt;
 547   size_t i;
 548
 549   if (lex_match_id ("ONEPASS"))
 550     two_pass = 0;
 551   if (token != '/')
 552     {
 553       lex_force_match ('/');
 554       goto done;
 555     }
 556   fprintf (stderr, "%s => ", lex_rest_of_line (NULL));
 557   lex_get ();
 558
 559   if (two_pass)
 560     {
 561       struct moments *m = NULL;
 562
 563       m = moments_create (MOMENT_KURTOSIS);
 564       if (!read_values (&values, &weights, &cnt))
 565         {
 566           moments_destroy (m);
 567           goto done;
 568         }
 569       for (i = 0; i < cnt; i++)
 570         moments_pass_one (m, values[i], weights[i]);
 571       for (i = 0; i < cnt; i++)
 572         moments_pass_two (m, values[i], weights[i]);
 573       moments_calculate (m, &weight, &M[0], &M[1], &M[2], &M[3]);
 574       moments_destroy (m);
 575     }
 576   else
 577     {
 578       struct moments1 *m = NULL;
 579
 580       m = moments1_create (MOMENT_KURTOSIS);
 581       if (!read_values (&values, &weights, &cnt))
 582         {
 583           moments1_destroy (m);
 584           goto done;
 585         }
 586       for (i = 0; i < cnt; i++)
 587         moments1_add (m, values[i], weights[i]);
 588       moments1_calculate (m, &weight, &M[0], &M[1], &M[2], &M[3]);
 589       moments1_destroy (m);
 590     }
 591
 592   fprintf (stderr, "W=%.3f", weight);
 593   for (i = 0; i < 4; i++)
 594     {
 595       fprintf (stderr, " M%d=", i + 1);
 596       if (M[i] == SYSMIS)
 597         fprintf (stderr, "sysmis");
 598       else if (fabs (M[i]) <= 0.0005)
 599         fprintf (stderr, "0.000");
 600       else
 601         fprintf (stderr, "%.3f", M[i]);
 602     }
 603   fprintf (stderr, "\n");
 604
 605   retval = lex_end_of_command ();
 606
 607  done:
 608   free (values);
 609   free (weights);
 610   return retval;
 611 }