Fix missing @clicksequence problem with older Texinfo versions.
[pspp-builds.git] / 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 }