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/value.h>
28 #define _(msgid) gettext (msgid)
30 /* Calculates variance, skewness, and kurtosis into *VARIANCE,
31 *SKEWNESS, and *KURTOSIS if they are non-null and not greater
32 moments than MAX_MOMENT. Accepts W as the total weight, D1 as
33 the total deviation from the estimated mean, and D2, D3, and
34 D4 as the sum of the squares, cubes, and 4th powers,
35 respectively, of the deviation from the estimated mean. */
37 calc_moments (enum moment max_moment,
38 double w, double d1, double d2, double d3, double d4,
39 double *variance, double *skewness, double *kurtosis)
43 if (max_moment >= MOMENT_VARIANCE && w > 1.)
45 double s2 = (d2 - pow2 (d1) / w) / (w - 1.);
49 /* From _SPSS Statistical Algorithms, 2nd ed.,
50 0-918469-89-9, section "DESCRIPTIVES". */
51 if (fabs (*variance) >= 1e-20)
53 if (max_moment >= MOMENT_SKEWNESS && skewness != NULL && w > 2.)
55 double s3 = s2 * sqrt (s2);
56 double g1 = (w * d3) / ((w - 1.0) * (w - 2.0) * s3);
60 if (max_moment >= MOMENT_KURTOSIS && kurtosis != NULL && w > 3.)
62 double den = (w - 2.) * (w - 3.) * pow2 (s2);
63 double g2 = (w * (w + 1) * d4 / (w - 1.) / den
64 - 3. * pow2 (d2) / den);
72 /* Two-pass moments. */
74 /* A set of two-pass moments. */
77 enum moment max_moment; /* Highest-order moment we're computing. */
78 int pass; /* Current pass (1 or 2). */
81 double w1; /* Total weight for pass 1, so far. */
82 double sum; /* Sum of values so far. */
83 double mean; /* Mean = sum / w1. */
86 double w2; /* Total weight for pass 2, so far. */
87 double d1; /* Sum of deviations from the mean. */
88 double d2; /* Sum of squared deviations from the mean. */
89 double d3; /* Sum of cubed deviations from the mean. */
90 double d4; /* Sum of (deviations from the mean)**4. */
93 /* Initializes moments M for calculating moment MAX_MOMENT and
96 init_moments (struct moments *m, enum moment max_moment)
99 assert (max_moment == MOMENT_MEAN || max_moment == MOMENT_VARIANCE
100 || max_moment == MOMENT_SKEWNESS || max_moment == MOMENT_KURTOSIS);
101 m->max_moment = max_moment;
105 /* Clears out a set of moments so that it can be reused for a new
106 set of values. The moments to be calculated are not changed. */
108 moments_clear (struct moments *m)
115 /* Creates and returns a data structure for calculating moment
116 MAX_MOMENT and lower moments on a data series. The user
117 should call moments_pass_one() for each value in the series,
118 then call moments_pass_two() for the same set of values in the
119 same order, then call moments_calculate() to obtain the
120 moments. The user may ask for the mean at any time during the
121 first pass (using moments_calculate()), but otherwise no
122 statistics may be requested until the end of the second pass.
123 Call moments_destroy() when the moments are no longer
126 moments_create (enum moment max_moment)
128 struct moments *m = xmalloc (sizeof *m);
129 init_moments (m, max_moment);
133 /* Adds VALUE with the given WEIGHT to the calculation of
134 moments for the first pass. */
136 moments_pass_one (struct moments *m, double value, double weight)
139 assert (m->pass == 1);
141 if (value != SYSMIS && weight > 0.)
143 m->sum += value * weight;
148 /* Adds VALUE with the given WEIGHT to the calculation of
149 moments for the second pass. */
151 moments_pass_two (struct moments *m, double value, double weight)
158 m->mean = m->w1 != 0. ? m->sum / m->w1 : 0.;
159 m->d1 = m->d2 = m->d3 = m->d4 = 0.;
162 if (value != SYSMIS && weight >= 0.)
164 double d = value - m->mean;
165 double d1_delta = d * weight;
167 if (m->max_moment >= MOMENT_VARIANCE)
169 double d2_delta = d1_delta * d;
171 if (m->max_moment >= MOMENT_SKEWNESS)
173 double d3_delta = d2_delta * d;
175 if (m->max_moment >= MOMENT_KURTOSIS)
177 double d4_delta = d3_delta * d;
186 /* Calculates moments based on the input data. Stores the total
187 weight in *WEIGHT, the mean in *MEAN, the variance in
188 *VARIANCE, the skewness in *SKEWNESS, and the kurtosis in
189 *KURTOSIS. Any of these result parameters may be null
190 pointers, in which case the values are not calculated. If any
191 result cannot be calculated, either because they are undefined
192 based on the input data or because their moments are higher
193 than the maximum requested on moments_create(), then SYSMIS is
194 stored into that result. */
196 moments_calculate (const struct moments *m,
198 double *mean, double *variance,
199 double *skewness, double *kurtosis)
205 if (variance != NULL)
207 if (skewness != NULL)
209 if (kurtosis != NULL)
215 /* How many passes so far? */
218 /* In the first pass we can only calculate the mean. */
219 if (mean != NULL && m->w1 > 0.)
220 *mean = m->sum / m->w1;
224 /* After the second pass we can calculate any stat. We
225 don't support "online" computation during the second
226 pass, so As a simple self-check, the total weight for
227 the passes must agree. */
228 assert (m->pass == 2);
229 assert (m->w1 == m->w2);
235 calc_moments (m->max_moment,
236 m->w2, m->d1, m->d2, m->d3, m->d4,
237 variance, skewness, kurtosis);
242 /* Destroys a set of moments. */
244 moments_destroy (struct moments *m)
249 /* Calculates the requested moments on the CNT values in ARRAY.
250 Each value is given a weight of 1. The total weight is stored
251 into *WEIGHT (trivially) and the mean, variance, skewness, and
252 kurtosis are stored into *MEAN, *VARIANCE, *SKEWNESS, and
253 *KURTOSIS, respectively. Any of the result pointers may be
254 null, in which case no value is stored. */
256 moments_of_doubles (const double *array, size_t cnt,
258 double *mean, double *variance,
259 double *skewness, double *kurtosis)
261 enum moment max_moment;
265 if (kurtosis != NULL)
266 max_moment = MOMENT_KURTOSIS;
267 else if (skewness != NULL)
268 max_moment = MOMENT_SKEWNESS;
269 else if (variance != NULL)
270 max_moment = MOMENT_VARIANCE;
272 max_moment = MOMENT_MEAN;
274 init_moments (&m, max_moment);
275 for (idx = 0; idx < cnt; idx++)
276 moments_pass_one (&m, array[idx], 1.);
277 for (idx = 0; idx < cnt; idx++)
278 moments_pass_two (&m, array[idx], 1.);
279 moments_calculate (&m, weight, mean, variance, skewness, kurtosis);
282 /* Calculates the requested moments on the CNT numeric values in
283 ARRAY. Each value is given a weight of 1. The total weight
284 is stored into *WEIGHT (trivially) and the mean, variance,
285 skewness, and kurtosis are stored into *MEAN, *VARIANCE,
286 *SKEWNESS, and *KURTOSIS, respectively. Any of the result
287 pointers may be null, in which case no value is stored. */
289 moments_of_values (const union value *array, size_t cnt,
291 double *mean, double *variance,
292 double *skewness, double *kurtosis)
294 enum moment max_moment;
298 if (kurtosis != NULL)
299 max_moment = MOMENT_KURTOSIS;
300 else if (skewness != NULL)
301 max_moment = MOMENT_SKEWNESS;
302 else if (variance != NULL)
303 max_moment = MOMENT_VARIANCE;
305 max_moment = MOMENT_MEAN;
307 init_moments (&m, max_moment);
308 for (idx = 0; idx < cnt; idx++)
309 moments_pass_one (&m, array[idx].f, 1.);
310 for (idx = 0; idx < cnt; idx++)
311 moments_pass_two (&m, array[idx].f, 1.);
312 moments_calculate (&m, weight, mean, variance, skewness, kurtosis);
315 /* One-pass moments. */
317 /* A set of one-pass moments. */
320 enum moment max_moment; /* Highest-order moment we're computing. */
321 double w; /* Total weight so far. */
322 double d1; /* Sum of deviations from the mean. */
323 double d2; /* Sum of squared deviations from the mean. */
324 double d3; /* Sum of cubed deviations from the mean. */
325 double d4; /* Sum of (deviations from the mean)**4. */
328 /* Initializes one-pass moments M for calculating moment
329 MAX_MOMENT and lower moments. */
331 init_moments1 (struct moments1 *m, enum moment max_moment)
334 assert (max_moment == MOMENT_MEAN || max_moment == MOMENT_VARIANCE
335 || max_moment == MOMENT_SKEWNESS || max_moment == MOMENT_KURTOSIS);
336 m->max_moment = max_moment;
340 /* Clears out a set of one-pass moments so that it can be reused
341 for a new set of values. The moments to be calculated are not
344 moments1_clear (struct moments1 *m)
347 m->d1 = m->d2 = m->d3 = m->d4 = 0.;
350 /* Creates and returns a data structure for calculating moment
351 MAX_MOMENT and lower moments on a data series in a single
352 pass. The user should call moments1_add() for each value in
353 the series. The user may call moments1_calculate() to obtain
354 the current moments at any time. Call moments1_destroy() when
355 the moments are no longer needed.
357 One-pass moments should only be used when two passes over the
358 data are impractical. */
360 moments1_create (enum moment max_moment)
362 struct moments1 *m = xmalloc (sizeof *m);
363 init_moments1 (m, max_moment);
367 /* Adds VALUE with the given WEIGHT to the calculation of
370 moments1_add (struct moments1 *m, double value, double weight)
374 if (value != SYSMIS && weight > 0.)
380 v1 = (weight / m->w) * (value - m->d1);
383 if (m->max_moment >= MOMENT_VARIANCE)
386 double w_prev_w = m->w * prev_w;
387 double prev_m2 = m->d2;
389 m->d2 += w_prev_w / weight * v2;
390 if (m->max_moment >= MOMENT_SKEWNESS)
392 double w2 = weight * weight;
394 double prev_m3 = m->d3;
396 m->d3 += (-3. * v1 * prev_m2
397 + w_prev_w / w2 * (m->w - 2. * weight) * v3);
398 if (m->max_moment >= MOMENT_KURTOSIS)
400 double w3 = w2 * weight;
403 m->d4 += (-4. * v1 * prev_m3
405 + ((pow2 (m->w) - 3. * weight * prev_w)
406 * v4 * w_prev_w / w3));
413 /* Calculates one-pass moments based on the input data. Stores
414 the total weight in *WEIGHT, the mean in *MEAN, the variance
415 in *VARIANCE, the skewness in *SKEWNESS, and the kurtosis in
416 *KURTOSIS. Any of these result parameters may be null
417 pointers, in which case the values are not calculated. If any
418 result cannot be calculated, either because they are undefined
419 based on the input data or because their moments are higher
420 than the maximum requested on moments_create(), then SYSMIS is
421 stored into that result. */
423 moments1_calculate (const struct moments1 *m,
425 double *mean, double *variance,
426 double *skewness, double *kurtosis)
432 if (variance != NULL)
434 if (skewness != NULL)
436 if (kurtosis != NULL)
447 calc_moments (m->max_moment,
448 m->w, 0., m->d2, m->d3, m->d4,
449 variance, skewness, kurtosis);
453 /* Destroy one-pass moments M. */
455 moments1_destroy (struct moments1 *m)
460 /* Returns the standard error of the skewness for the given total
463 From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
464 section "DESCRIPTIVES". */
466 calc_seskew (double W)
468 return sqrt ((6. * W * (W - 1.)) / ((W - 2.) * (W + 1.) * (W + 3.)));
471 /* Returns the standard error of the kurtosis for the given total
474 From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
475 section "DESCRIPTIVES", except that the sqrt symbol is omitted
478 calc_sekurt (double W)
480 return sqrt ((4. * (pow2 (W) - 1.) * pow2 (calc_seskew (W)))
481 / ((W - 3.) * (W + 5.)));