gui: Fix behavior of print preview in psppire-output-window.
[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.  We
229          don't support "online" computation during the second
230          pass, so As a simple self-check, the total weight for
231          the passes must agree. */
232       assert (m->pass == 2);
233       assert (m->w1 == m->w2);
234
235       if (m->w2 > 0.)
236         {
237           if (mean != NULL)
238             *mean = m->mean;
239           calc_moments (m->max_moment,
240                         m->w2, m->d1, m->d2, m->d3, m->d4,
241                         variance, skewness, kurtosis);
242         }
243     }
244 }
245
246 /* Destroys a set of moments. */
247 void
248 moments_destroy (struct moments *m)
249 {
250   free (m);
251 }
252
253 /* Calculates the requested moments on the CNT values in ARRAY.
254    Each value is given a weight of 1.  The total weight is stored
255    into *WEIGHT (trivially) and the mean, variance, skewness, and
256    kurtosis are stored into *MEAN, *VARIANCE, *SKEWNESS, and
257    *KURTOSIS, respectively.  Any of the result pointers may be
258    null, in which case no value is stored. */
259 void
260 moments_of_doubles (const double *array, size_t cnt,
261                     double *weight,
262                     double *mean, double *variance,
263                     double *skewness, double *kurtosis)
264 {
265   enum moment max_moment;
266   struct moments m;
267   size_t idx;
268
269   if (kurtosis != NULL)
270     max_moment = MOMENT_KURTOSIS;
271   else if (skewness != NULL)
272     max_moment = MOMENT_SKEWNESS;
273   else if (variance != NULL)
274     max_moment = MOMENT_VARIANCE;
275   else
276     max_moment = MOMENT_MEAN;
277
278   init_moments (&m, max_moment);
279   for (idx = 0; idx < cnt; idx++)
280     moments_pass_one (&m, array[idx], 1.);
281   for (idx = 0; idx < cnt; idx++)
282     moments_pass_two (&m, array[idx], 1.);
283   moments_calculate (&m, weight, mean, variance, skewness, kurtosis);
284 }
285
286 /* Calculates the requested moments on the CNT numeric values in
287    ARRAY.  Each value is given a weight of 1.  The total weight
288    is stored into *WEIGHT (trivially) and the mean, variance,
289    skewness, and kurtosis are stored into *MEAN, *VARIANCE,
290    *SKEWNESS, and *KURTOSIS, respectively.  Any of the result
291    pointers may be null, in which case no value is stored. */
292 void
293 moments_of_values (const union value *array, size_t cnt,
294                    double *weight,
295                    double *mean, double *variance,
296                    double *skewness, double *kurtosis)
297 {
298   enum moment max_moment;
299   struct moments m;
300   size_t idx;
301
302   if (kurtosis != NULL)
303     max_moment = MOMENT_KURTOSIS;
304   else if (skewness != NULL)
305     max_moment = MOMENT_SKEWNESS;
306   else if (variance != NULL)
307     max_moment = MOMENT_VARIANCE;
308   else
309     max_moment = MOMENT_MEAN;
310
311   init_moments (&m, max_moment);
312   for (idx = 0; idx < cnt; idx++)
313     moments_pass_one (&m, array[idx].f, 1.);
314   for (idx = 0; idx < cnt; idx++)
315     moments_pass_two (&m, array[idx].f, 1.);
316   moments_calculate (&m, weight, mean, variance, skewness, kurtosis);
317 }
318 \f
319 /* One-pass moments. */
320
321 /* A set of one-pass moments. */
322 struct moments1
323   {
324     enum moment max_moment;     /* Highest-order moment we're computing. */
325     double w;                   /* Total weight so far. */
326     double d1;                  /* Sum of deviations from the mean. */
327     double d2;                  /* Sum of squared deviations from the mean. */
328     double d3;                  /* Sum of cubed deviations from the mean. */
329     double d4;                  /* Sum of (deviations from the mean)**4. */
330   };
331
332 /* Initializes one-pass moments M for calculating moment
333    MAX_MOMENT and lower moments. */
334 static void
335 init_moments1 (struct moments1 *m, enum moment max_moment)
336 {
337   assert (m != NULL);
338   assert (max_moment == MOMENT_MEAN || max_moment == MOMENT_VARIANCE
339           || max_moment == MOMENT_SKEWNESS || max_moment == MOMENT_KURTOSIS);
340   m->max_moment = max_moment;
341   moments1_clear (m);
342 }
343
344 /* Clears out a set of one-pass moments so that it can be reused
345    for a new set of values.  The moments to be calculated are not
346    changed. */
347 void
348 moments1_clear (struct moments1 *m)
349 {
350   m->w = 0.;
351   m->d1 = m->d2 = m->d3 = m->d4 = 0.;
352 }
353
354 /* Creates and returns a data structure for calculating moment
355    MAX_MOMENT and lower moments on a data series in a single
356    pass.  The user should call moments1_add() for each value in
357    the series.  The user may call moments1_calculate() to obtain
358    the current moments at any time.  Call moments1_destroy() when
359    the moments are no longer needed.
360
361    One-pass moments should only be used when two passes over the
362    data are impractical. */
363 struct moments1 *
364 moments1_create (enum moment max_moment)
365 {
366   struct moments1 *m = xmalloc (sizeof *m);
367   init_moments1 (m, max_moment);
368   return m;
369 }
370
371 /* Adds VALUE with the given WEIGHT to the calculation of
372    one-pass moments. */
373 void
374 moments1_add (struct moments1 *m, double value, double weight)
375 {
376   assert (m != NULL);
377
378   if (value != SYSMIS && weight > 0.)
379     {
380       double prev_w, v1;
381
382       prev_w = m->w;
383       m->w += weight;
384       v1 = (weight / m->w) * (value - m->d1);
385       m->d1 += v1;
386
387       if (m->max_moment >= MOMENT_VARIANCE)
388         {
389           double v2 = v1 * v1;
390           double w_prev_w = m->w * prev_w;
391           double prev_m2 = m->d2;
392
393           m->d2 += w_prev_w / weight * v2;
394           if (m->max_moment >= MOMENT_SKEWNESS)
395             {
396               double w2 = weight * weight;
397               double v3 = v2 * v1;
398               double prev_m3 = m->d3;
399
400               m->d3 += (-3. * v1 * prev_m2
401                          + w_prev_w / w2 * (m->w - 2. * weight) * v3);
402               if (m->max_moment >= MOMENT_KURTOSIS)
403                 {
404                   double w3 = w2 * weight;
405                   double v4 = v2 * v2;
406
407                   m->d4 += (-4. * v1 * prev_m3
408                              + 6. * v2 * prev_m2
409                              + ((pow2 (m->w) - 3. * weight * prev_w)
410                                 * v4 * w_prev_w / w3));
411                 }
412             }
413         }
414     }
415 }
416
417 /* Calculates one-pass moments based on the input data.  Stores
418    the total weight in *WEIGHT, the mean in *MEAN, the variance
419    in *VARIANCE, the skewness in *SKEWNESS, and the kurtosis in
420    *KURTOSIS.  Any of these result parameters may be null
421    pointers, in which case the values are not calculated.  If any
422    result cannot be calculated, either because they are undefined
423    based on the input data or because their moments are higher
424    than the maximum requested on moments_create(), then SYSMIS is
425    stored into that result. */
426 void
427 moments1_calculate (const struct moments1 *m,
428                     double *weight,
429                     double *mean, double *variance,
430                     double *skewness, double *kurtosis)
431 {
432   assert (m != NULL);
433
434   if (mean != NULL)
435     *mean = SYSMIS;
436   if (variance != NULL)
437     *variance = SYSMIS;
438   if (skewness != NULL)
439     *skewness = SYSMIS;
440   if (kurtosis != NULL)
441     *kurtosis = SYSMIS;
442
443   if (weight != NULL)
444     *weight = m->w;
445
446   if (m->w > 0.)
447     {
448       if (mean != NULL)
449         *mean = m->d1;
450
451       calc_moments (m->max_moment,
452                     m->w, 0., m->d2, m->d3, m->d4,
453                     variance, skewness, kurtosis);
454     }
455 }
456
457 /* Destroy one-pass moments M. */
458 void
459 moments1_destroy (struct moments1 *m)
460 {
461   free (m);
462 }
463 \f
464 /* Returns the standard error of the skewness for the given total
465    weight W.
466
467    From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
468    section "DESCRIPTIVES". */
469 double
470 calc_seskew (double W)
471 {
472   return sqrt ((6. * W * (W - 1.)) / ((W - 2.) * (W + 1.) * (W + 3.)));
473 }
474
475 /* Returns the standard error of the kurtosis for the given total
476    weight W.
477
478    From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
479    section "DESCRIPTIVES", except that the sqrt symbol is omitted
480    there. */
481 double
482 calc_sekurt (double W)
483 {
484   return sqrt ((4. * (pow2 (W) - 1.) * pow2 (calc_seskew (W)))
485                / ((W - 3.) * (W + 5.)));
486 }