src/math/moments.c

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