src/math/moments.c

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