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