gui: Allow File|Open to select an encoding for system files.
[pspp] / 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.  */
229       assert (m->pass == 2);
230
231       if (m->w2 > 0.)
232         {
233           if (mean != NULL)
234             *mean = m->mean;
235           calc_moments (m->max_moment,
236                         m->w2, m->d1, m->d2, m->d3, m->d4,
237                         variance, skewness, kurtosis);
238         }
239     }
240 }
241
242 /* Destroys a set of moments. */
243 void
244 moments_destroy (struct moments *m)
245 {
246   free (m);
247 }
248
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. */
255 void
256 moments_of_doubles (const double *array, size_t cnt,
257                     double *weight,
258                     double *mean, double *variance,
259                     double *skewness, double *kurtosis)
260 {
261   enum moment max_moment;
262   struct moments m;
263   size_t idx;
264
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;
271   else
272     max_moment = MOMENT_MEAN;
273
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);
280 }
281
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. */
288 void
289 moments_of_values (const union value *array, size_t cnt,
290                    double *weight,
291                    double *mean, double *variance,
292                    double *skewness, double *kurtosis)
293 {
294   enum moment max_moment;
295   struct moments m;
296   size_t idx;
297
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;
304   else
305     max_moment = MOMENT_MEAN;
306
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);
313 }
314 \f
315 /* One-pass moments. */
316
317 /* A set of one-pass moments. */
318 struct moments1
319   {
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. */
326   };
327
328 /* Initializes one-pass moments M for calculating moment
329    MAX_MOMENT and lower moments. */
330 static void
331 init_moments1 (struct moments1 *m, enum moment max_moment)
332 {
333   assert (m != NULL);
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;
337   moments1_clear (m);
338 }
339
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
342    changed. */
343 void
344 moments1_clear (struct moments1 *m)
345 {
346   m->w = 0.;
347   m->d1 = m->d2 = m->d3 = m->d4 = 0.;
348 }
349
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.
356
357    One-pass moments should only be used when two passes over the
358    data are impractical. */
359 struct moments1 *
360 moments1_create (enum moment max_moment)
361 {
362   struct moments1 *m = xmalloc (sizeof *m);
363   init_moments1 (m, max_moment);
364   return m;
365 }
366
367 /* Adds VALUE with the given WEIGHT to the calculation of
368    one-pass moments. */
369 void
370 moments1_add (struct moments1 *m, double value, double weight)
371 {
372   assert (m != NULL);
373
374   if (value != SYSMIS && weight > 0.)
375     {
376       double prev_w, v1;
377
378       prev_w = m->w;
379       m->w += weight;
380       v1 = (weight / m->w) * (value - m->d1);
381       m->d1 += v1;
382
383       if (m->max_moment >= MOMENT_VARIANCE)
384         {
385           double v2 = v1 * v1;
386           double w_prev_w = m->w * prev_w;
387           double prev_m2 = m->d2;
388
389           m->d2 += w_prev_w / weight * v2;
390           if (m->max_moment >= MOMENT_SKEWNESS)
391             {
392               double w2 = weight * weight;
393               double v3 = v2 * v1;
394               double prev_m3 = m->d3;
395
396               m->d3 += (-3. * v1 * prev_m2
397                          + w_prev_w / w2 * (m->w - 2. * weight) * v3);
398               if (m->max_moment >= MOMENT_KURTOSIS)
399                 {
400                   double w3 = w2 * weight;
401                   double v4 = v2 * v2;
402
403                   m->d4 += (-4. * v1 * prev_m3
404                              + 6. * v2 * prev_m2
405                              + ((pow2 (m->w) - 3. * weight * prev_w)
406                                 * v4 * w_prev_w / w3));
407                 }
408             }
409         }
410     }
411 }
412
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. */
422 void
423 moments1_calculate (const struct moments1 *m,
424                     double *weight,
425                     double *mean, double *variance,
426                     double *skewness, double *kurtosis)
427 {
428   assert (m != NULL);
429
430   if (mean != NULL)
431     *mean = SYSMIS;
432   if (variance != NULL)
433     *variance = SYSMIS;
434   if (skewness != NULL)
435     *skewness = SYSMIS;
436   if (kurtosis != NULL)
437     *kurtosis = SYSMIS;
438
439   if (weight != NULL)
440     *weight = m->w;
441
442   if (m->w > 0.)
443     {
444       if (mean != NULL)
445         *mean = m->d1;
446
447       calc_moments (m->max_moment,
448                     m->w, 0., m->d2, m->d3, m->d4,
449                     variance, skewness, kurtosis);
450     }
451 }
452
453 /* Destroy one-pass moments M. */
454 void
455 moments1_destroy (struct moments1 *m)
456 {
457   free (m);
458 }
459 \f
460
461 double
462 calc_semean (double var, double W)
463 {
464   return sqrt (var / W);
465 }
466
467 /* Returns the standard error of the skewness for the given total
468    weight W.
469
470    From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
471    section "DESCRIPTIVES". */
472 double
473 calc_seskew (double W)
474 {
475   return sqrt ((6. * W * (W - 1.)) / ((W - 2.) * (W + 1.) * (W + 3.)));
476 }
477
478 /* Returns the standard error of the kurtosis for the given total
479    weight W.
480
481    From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
482    section "DESCRIPTIVES", except that the sqrt symbol is omitted
483    there. */
484 double
485 calc_sekurt (double W)
486 {
487   return sqrt ((4. * (pow2 (W) - 1.) * pow2 (calc_seskew (W)))
488                / ((W - 3.) * (W + 5.)));
489 }