1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
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.
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.
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/>. */
22 #include <libpspp/misc.h>
23 #include <data/val-type.h>
24 #include <data/value.h>
29 #define _(msgid) gettext (msgid)
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. */
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)
44 if (max_moment >= MOMENT_VARIANCE && w > 1.)
46 double s2 = (d2 - pow2 (d1) / w) / (w - 1.);
50 /* From _SPSS Statistical Algorithms, 2nd ed.,
51 0-918469-89-9, section "DESCRIPTIVES". */
52 if (fabs (*variance) >= 1e-20)
54 if (max_moment >= MOMENT_SKEWNESS && skewness != NULL && w > 2.)
56 double s3 = s2 * sqrt (s2);
57 double g1 = (w * d3) / ((w - 1.0) * (w - 2.0) * s3);
61 if (max_moment >= MOMENT_KURTOSIS && kurtosis != NULL && w > 3.)
63 double den = (w - 2.) * (w - 3.) * pow2 (s2);
64 double g2 = (w * (w + 1) * d4 / (w - 1.) / den
65 - 3. * pow2 (d2) / den);
73 /* Two-pass moments. */
75 /* A set of two-pass moments. */
78 enum moment max_moment; /* Highest-order moment we're computing. */
79 int pass; /* Current pass (1 or 2). */
82 double w1; /* Total weight for pass 1, so far. */
83 double sum; /* Sum of values so far. */
84 double mean; /* Mean = sum / w1. */
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. */
94 /* Initializes moments M for calculating moment MAX_MOMENT and
97 init_moments (struct moments *m, enum moment max_moment)
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;
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. */
109 moments_clear (struct moments *m)
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
127 moments_create (enum moment max_moment)
129 struct moments *m = xmalloc (sizeof *m);
130 init_moments (m, max_moment);
134 /* Adds VALUE with the given WEIGHT to the calculation of
135 moments for the first pass. */
137 moments_pass_one (struct moments *m, double value, double weight)
140 assert (m->pass == 1);
142 if (value != SYSMIS && weight > 0.)
144 m->sum += value * weight;
149 /* Adds VALUE with the given WEIGHT to the calculation of
150 moments for the second pass. */
152 moments_pass_two (struct moments *m, double value, double weight)
159 m->mean = m->w1 != 0. ? m->sum / m->w1 : 0.;
160 m->d1 = m->d2 = m->d3 = m->d4 = 0.;
163 if (value != SYSMIS && weight >= 0.)
165 double d = value - m->mean;
166 double d1_delta = d * weight;
168 if (m->max_moment >= MOMENT_VARIANCE)
170 double d2_delta = d1_delta * d;
172 if (m->max_moment >= MOMENT_SKEWNESS)
174 double d3_delta = d2_delta * d;
176 if (m->max_moment >= MOMENT_KURTOSIS)
178 double d4_delta = d3_delta * d;
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. */
197 moments_calculate (const struct moments *m,
199 double *mean, double *variance,
200 double *skewness, double *kurtosis)
206 if (variance != NULL)
208 if (skewness != NULL)
210 if (kurtosis != NULL)
216 /* How many passes so far? */
219 /* In the first pass we can only calculate the mean. */
220 if (mean != NULL && m->w1 > 0.)
221 *mean = m->sum / m->w1;
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);
236 calc_moments (m->max_moment,
237 m->w2, m->d1, m->d2, m->d3, m->d4,
238 variance, skewness, kurtosis);
243 /* Destroys a set of moments. */
245 moments_destroy (struct moments *m)
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. */
257 moments_of_doubles (const double *array, size_t cnt,
259 double *mean, double *variance,
260 double *skewness, double *kurtosis)
262 enum moment max_moment;
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;
273 max_moment = MOMENT_MEAN;
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);
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. */
290 moments_of_values (const union value *array, size_t cnt,
292 double *mean, double *variance,
293 double *skewness, double *kurtosis)
295 enum moment max_moment;
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;
306 max_moment = MOMENT_MEAN;
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);
316 /* One-pass moments. */
318 /* A set of one-pass moments. */
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. */
329 /* Initializes one-pass moments M for calculating moment
330 MAX_MOMENT and lower moments. */
332 init_moments1 (struct moments1 *m, enum moment max_moment)
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;
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
345 moments1_clear (struct moments1 *m)
348 m->d1 = m->d2 = m->d3 = m->d4 = 0.;
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.
358 One-pass moments should only be used when two passes over the
359 data are impractical. */
361 moments1_create (enum moment max_moment)
363 struct moments1 *m = xmalloc (sizeof *m);
364 init_moments1 (m, max_moment);
368 /* Adds VALUE with the given WEIGHT to the calculation of
371 moments1_add (struct moments1 *m, double value, double weight)
375 if (value != SYSMIS && weight > 0.)
381 v1 = (weight / m->w) * (value - m->d1);
384 if (m->max_moment >= MOMENT_VARIANCE)
387 double w_prev_w = m->w * prev_w;
388 double prev_m2 = m->d2;
390 m->d2 += w_prev_w / weight * v2;
391 if (m->max_moment >= MOMENT_SKEWNESS)
393 double w2 = weight * weight;
395 double prev_m3 = m->d3;
397 m->d3 += (-3. * v1 * prev_m2
398 + w_prev_w / w2 * (m->w - 2. * weight) * v3);
399 if (m->max_moment >= MOMENT_KURTOSIS)
401 double w3 = w2 * weight;
404 m->d4 += (-4. * v1 * prev_m3
406 + ((pow2 (m->w) - 3. * weight * prev_w)
407 * v4 * w_prev_w / w3));
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. */
424 moments1_calculate (const struct moments1 *m,
426 double *mean, double *variance,
427 double *skewness, double *kurtosis)
433 if (variance != NULL)
435 if (skewness != NULL)
437 if (kurtosis != NULL)
448 calc_moments (m->max_moment,
449 m->w, 0., m->d2, m->d3, m->d4,
450 variance, skewness, kurtosis);
454 /* Destroy one-pass moments M. */
456 moments1_destroy (struct moments1 *m)
461 /* Returns the standard error of the skewness for the given total
464 From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
465 section "DESCRIPTIVES". */
467 calc_seskew (double W)
469 return sqrt ((6. * W * (W - 1.)) / ((W - 2.) * (W + 1.) * (W + 3.)));
472 /* Returns the standard error of the kurtosis for the given total
475 From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
476 section "DESCRIPTIVES", except that the sqrt symbol is omitted
479 calc_sekurt (double W)
481 return sqrt ((4. * (pow2 (W) - 1.) * pow2 (calc_seskew (W)))
482 / ((W - 3.) * (W + 5.)));