Fix assertion for proper Huffman merge pattern: 0 == 1 modulo 1.
[pspp] / src / moments.c
1 /* PSPP - computes sample statistics.
2    Copyright (C) 1997-9, 2000 Free Software Foundation, Inc.
3    Written by Ben Pfaff <blp@gnu.org>.
4
5    This program is free software; you can redistribute it and/or
6    modify it under the terms of the GNU General Public License as
7    published by the Free Software Foundation; either version 2 of the
8    License, or (at your option) any later version.
9
10    This program is distributed in the hope that it will be useful, but
11    WITHOUT ANY WARRANTY; without even the implied warranty of
12    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
13    General Public License for more details.
14
15    You should have received a copy of the GNU General Public License
16    along with this program; if not, write to the Free Software
17    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
18    02111-1307, USA. */
19
20 #include <config.h>
21 #include "moments.h"
22 #include <assert.h>
23 #include <math.h>
24 #include <stdlib.h>
25 #include "alloc.h"
26 #include "misc.h"
27 #include "val.h"
28 \f
29 /* Calculates variance, skewness, and kurtosis into *VARIANCE,
30    *SKEWNESS, and *KURTOSIS if they are non-null and not greater
31    moments than MAX_MOMENT.  Accepts W as the total weight, D1 as
32    the total deviation from the estimated mean, and D2, D3, and
33    D4 as the sum of the squares, cubes, and 4th powers,
34    respectively, of the deviation from the estimated mean. */
35 static void
36 calc_moments (enum moment max_moment,
37               double w, double d1, double d2, double d3, double d4,
38               double *variance, double *skewness, double *kurtosis) 
39 {
40   assert (w > 0.);
41
42   if (max_moment >= MOMENT_VARIANCE && w > 1.) 
43     {
44       double s2;
45
46       /* From _Numerical Recipes in C_, 2nd ed., 0-521-43108-5,
47          section 14.1. */
48       s2 = (d2 - pow2 (d1) / w) / (w - 1.);
49       if (variance != NULL)
50         *variance = s2;
51
52       /* From _SPSS Statistical Algorithms, 2nd ed.,
53          0-918469-89-9, section "DESCRIPTIVES". */
54       if (fabs (*variance) >= 1e-20) 
55         {
56           if (max_moment >= MOMENT_SKEWNESS && skewness != NULL && w > 2.)
57             {
58               double s3 = s2 * sqrt (s2);
59               double g1 = (w * d3) / ((w - 1.0) * (w - 2.0) * s3);
60               if (finite (g1))
61                 *skewness = g1; 
62             }
63           if (max_moment >= MOMENT_KURTOSIS && kurtosis != NULL && w > 3.)
64             {
65               double den = (w - 2.) * (w - 3.) * pow2 (s2);
66               double g2 = (w * (w + 1) * d4 / (w - 1.) / den
67                            - 3. * pow2 (d2) / den);
68               if (finite (g2))
69                 *kurtosis = g2; 
70             }
71         } 
72     }
73 }
74 \f
75 /* Two-pass moments. */
76
77 /* A set of two-pass moments. */
78 struct moments 
79   {
80     enum moment max_moment;     /* Highest-order moment we're computing. */
81     int pass;                   /* Current pass (1 or 2). */
82
83     /* Pass one. */
84     double w1;                  /* Total weight for pass 1, so far. */
85     double sum;                 /* Sum of values so far. */
86     double mean;                /* Mean = sum / w1. */
87
88     /* Pass two. */
89     double w2;                  /* Total weight for pass 2, so far. */
90     double d1;                  /* Sum of deviations from the mean. */
91     double d2;                  /* Sum of squared deviations from the mean. */
92     double d3;                  /* Sum of cubed deviations from the mean. */
93     double d4;                  /* Sum of (deviations from the mean)**4. */
94   };
95
96 /* Initializes moments M for calculating moment MAX_MOMENT and
97    lower moments. */
98 static void
99 init_moments (struct moments *m, enum moment max_moment)
100 {
101   assert (m != NULL);
102   assert (max_moment == MOMENT_MEAN || max_moment == MOMENT_VARIANCE
103           || max_moment == MOMENT_SKEWNESS || max_moment == MOMENT_KURTOSIS);
104   m->max_moment = max_moment;
105   moments_clear (m);
106 }
107
108 /* Clears out a set of moments so that it can be reused for a new
109    set of values.  The moments to be calculated are not changed. */
110 void
111 moments_clear (struct moments *m) 
112 {
113   m->pass = 1;
114   m->w1 = m->w2 = 0.;
115   m->sum = 0.;
116 }
117
118 /* Creates and returns a data structure for calculating moment
119    MAX_MOMENT and lower moments on a data series.  The user
120    should call moments_pass_one() for each value in the series,
121    then call moments_pass_two() for the same set of values in the
122    same order, then call moments_calculate() to obtain the
123    moments.  The user may ask for the mean at any time during the
124    first pass (using moments_calculate()), but otherwise no
125    statistics may be requested until the end of the second pass.
126    Call moments_destroy() when the moments are no longer
127    needed. */
128 struct moments *
129 moments_create (enum moment max_moment)
130 {
131   struct moments *m = xmalloc (sizeof *m);
132   init_moments (m, max_moment);
133   return m;
134 }
135
136 /* Adds VALUE with the given WEIGHT to the calculation of
137    moments for the first pass. */
138 void
139 moments_pass_one (struct moments *m, double value, double weight) 
140 {
141   assert (m != NULL);
142   assert (m->pass == 1);
143
144   if (value != SYSMIS && weight > 0.) 
145     {
146       m->sum += value * weight;
147       m->w1 += weight;
148     }
149 }
150
151 /* Adds VALUE with the given WEIGHT to the calculation of
152    moments for the second pass. */
153 void
154 moments_pass_two (struct moments *m, double value, double weight) 
155 {
156   double d, d_power;
157
158   assert (m != NULL);
159
160   if (m->pass == 1) 
161     {
162       m->pass = 2;
163       m->mean = m->w1 != 0. ? m->sum / m->w1 : 0.;
164       m->d1 = m->d2 = m->d3 = m->d4 = 0.;
165     }
166
167   if (value != SYSMIS && weight >= 0.) 
168     {
169       m->w2 += weight;
170
171       d = d_power = value - m->mean;
172       m->d1 += d_power * weight;
173
174       if (m->max_moment >= MOMENT_VARIANCE) 
175         {
176           d_power *= d;
177           m->d2 += d_power * weight;
178
179           if (m->max_moment >= MOMENT_SKEWNESS)
180             {
181               d_power *= d;
182               m->d3 += d_power * weight;
183
184               if (m->max_moment >= MOMENT_KURTOSIS)
185                 {
186                   d_power *= d;
187                   m->d4 += d_power * weight;
188                 }
189             }
190         }
191     }
192 }
193
194 /* Calculates moments based on the input data.  Stores the total
195    weight in *WEIGHT, the mean in *MEAN, the variance in
196    *VARIANCE, the skewness in *SKEWNESS, and the kurtosis in
197    *KURTOSIS.  Any of these result parameters may be null
198    pointers, in which case the values are not calculated.  If any
199    result cannot be calculated, either because they are undefined
200    based on the input data or because their moments are higher
201    than the maximum requested on moments_create(), then SYSMIS is
202    stored into that result. */
203 void
204 moments_calculate (const struct moments *m,
205                    double *weight,
206                    double *mean, double *variance,
207                    double *skewness, double *kurtosis) 
208 {
209   assert (m != NULL);
210
211   if (mean != NULL)
212     *mean = SYSMIS;
213   if (variance != NULL)
214     *variance = SYSMIS;
215   if (skewness != NULL)
216     *skewness = SYSMIS;
217   if (kurtosis != NULL)
218     *kurtosis = SYSMIS;
219
220   if (weight != NULL)
221     *weight = m->w1;
222
223   /* How many passes so far? */
224   if (m->pass == 1) 
225     {
226       /* In the first pass we can only calculate the mean. */
227       if (mean != NULL && m->w1 > 0.)
228         *mean = m->sum / m->w1;
229     }
230   else 
231     {
232       /* After the second pass we can calculate any stat.  We
233          don't support "online" computation during the second
234          pass, so As a simple self-check, the total weight for
235          the passes must agree. */
236       assert (m->pass == 2);
237       assert (m->w1 == m->w2);
238
239       if (m->w2 > 0.) 
240         {
241           if (mean != NULL)
242             *mean = m->mean;
243           calc_moments (m->max_moment,
244                         m->w2, m->d1, m->d2, m->d3, m->d4,
245                         variance, skewness, kurtosis); 
246         }
247     }
248 }
249
250 /* Destroys a set of moments. */
251 void
252 moments_destroy (struct moments *m) 
253 {
254   free (m);
255 }
256
257 /* Calculates the requested moments on the CNT values in ARRAY.
258    Each value is given a weight of 1.  The total weight is stored
259    into *WEIGHT (trivially) and the mean, variance, skewness, and
260    kurtosis are stored into *MEAN, *VARIANCE, *SKEWNESS, and
261    *KURTOSIS, respectively.  Any of the result pointers may be
262    null, in which case no value is stored. */
263 void
264 moments_of_doubles (const double *array, size_t cnt,
265                     double *weight,
266                     double *mean, double *variance,
267                     double *skewness, double *kurtosis) 
268 {
269   enum moment max_moment;
270   struct moments m;
271   size_t idx;
272
273   if (kurtosis != NULL)
274     max_moment = MOMENT_KURTOSIS;
275   else if (skewness != NULL)
276     max_moment = MOMENT_SKEWNESS;
277   else if (variance != NULL)
278     max_moment = MOMENT_VARIANCE;
279   else
280     max_moment = MOMENT_MEAN;
281
282   init_moments (&m, max_moment);
283   for (idx = 0; idx < cnt; idx++)
284     moments_pass_one (&m, array[idx], 1.);
285   for (idx = 0; idx < cnt; idx++)
286     moments_pass_two (&m, array[idx], 1.);
287   moments_calculate (&m, weight, mean, variance, skewness, kurtosis);
288 }
289
290 /* Calculates the requested moments on the CNT numeric values in
291    ARRAY.  Each value is given a weight of 1.  The total weight
292    is stored into *WEIGHT (trivially) and the mean, variance,
293    skewness, and kurtosis are stored into *MEAN, *VARIANCE,
294    *SKEWNESS, and *KURTOSIS, respectively.  Any of the result
295    pointers may be null, in which case no value is stored. */
296 void
297 moments_of_values (const union value *array, size_t cnt,
298                    double *weight,
299                    double *mean, double *variance,
300                    double *skewness, double *kurtosis) 
301 {
302   enum moment max_moment;
303   struct moments m;
304   size_t idx;
305
306   if (kurtosis != NULL)
307     max_moment = MOMENT_KURTOSIS;
308   else if (skewness != NULL)
309     max_moment = MOMENT_SKEWNESS;
310   else if (variance != NULL)
311     max_moment = MOMENT_VARIANCE;
312   else
313     max_moment = MOMENT_MEAN;
314
315   init_moments (&m, max_moment);
316   for (idx = 0; idx < cnt; idx++)
317     moments_pass_one (&m, array[idx].f, 1.);
318   for (idx = 0; idx < cnt; idx++)
319     moments_pass_two (&m, array[idx].f, 1.);
320   moments_calculate (&m, weight, mean, variance, skewness, kurtosis);
321 }
322 \f
323 /* One-pass moments. */
324
325 /* A set of one-pass moments. */
326 struct moments1 
327   {
328     enum moment max_moment;     /* Highest-order moment we're computing. */
329     double w;                   /* Total weight so far. */
330     double d1;                  /* Sum of deviations from the mean. */
331     double d2;                  /* Sum of squared deviations from the mean. */
332     double d3;                  /* Sum of cubed deviations from the mean. */
333     double d4;                  /* Sum of (deviations from the mean)**4. */
334   };
335
336 /* Initializes one-pass moments M for calculating moment
337    MAX_MOMENT and lower moments. */
338 static void
339 init_moments1 (struct moments1 *m, enum moment max_moment)
340 {
341   assert (m != NULL);
342   assert (max_moment == MOMENT_MEAN || max_moment == MOMENT_VARIANCE
343           || max_moment == MOMENT_SKEWNESS || max_moment == MOMENT_KURTOSIS);
344   m->max_moment = max_moment;
345   moments1_clear (m);
346 }
347
348 /* Clears out a set of one-pass moments so that it can be reused
349    for a new set of values.  The moments to be calculated are not
350    changed. */
351 void
352 moments1_clear (struct moments1 *m) 
353 {
354   m->w = 0.;
355   m->d1 = m->d2 = m->d3 = m->d4 = 0.;
356 }
357
358 /* Creates and returns a data structure for calculating moment
359    MAX_MOMENT and lower moments on a data series in a single
360    pass.  The user should call moments1_add() for each value in
361    the series.  The user may call moments1_calculate() to obtain
362    the current moments at any time.  Call moments1_destroy() when
363    the moments are no longer needed. 
364
365    One-pass moments should only be used when two passes over the
366    data are impractical. */
367 struct moments1 *
368 moments1_create (enum moment max_moment) 
369 {
370   struct moments1 *m = xmalloc (sizeof *m);
371   init_moments1 (m, max_moment);
372   return m;
373 }
374
375 /* Adds VALUE with the given WEIGHT to the calculation of
376    one-pass moments. */
377 void
378 moments1_add (struct moments1 *m, double value, double weight) 
379 {
380   assert (m != NULL);
381
382   if (value != SYSMIS && weight > 0.) 
383     {
384       double prev_w, v1;
385
386       prev_w = m->w;
387       m->w += weight;
388       v1 = (weight / m->w) * (value - m->d1);
389       m->d1 += v1;
390
391       if (m->max_moment >= MOMENT_VARIANCE) 
392         {
393           double v2 = v1 * v1;
394           double w_prev_w = m->w * prev_w;
395           double prev_m2 = m->d2;
396           
397           m->d2 += w_prev_w / weight * v2;
398           if (m->max_moment >= MOMENT_SKEWNESS) 
399             {
400               double w2 = weight * weight;
401               double v3 = v2 * v1;
402               double prev_m3 = m->d3;
403
404               m->d3 += (-3. * v1 * prev_m2
405                          + w_prev_w / w2 * (m->w - 2. * weight) * v3);
406               if (m->max_moment >= MOMENT_KURTOSIS) 
407                 {
408                   double w3 = w2 * weight;
409                   double v4 = v2 * v2;
410
411                   m->d4 += (-4. * v1 * prev_m3
412                              + 6. * v2 * prev_m2
413                              + ((pow2 (m->w) - 3. * weight * prev_w)
414                                 * v4 * w_prev_w / w3));
415                 }
416             }
417         }
418     }
419 }
420
421 /* Calculates one-pass moments based on the input data.  Stores
422    the total weight in *WEIGHT, the mean in *MEAN, the variance
423    in *VARIANCE, the skewness in *SKEWNESS, and the kurtosis in
424    *KURTOSIS.  Any of these result parameters may be null
425    pointers, in which case the values are not calculated.  If any
426    result cannot be calculated, either because they are undefined
427    based on the input data or because their moments are higher
428    than the maximum requested on moments_create(), then SYSMIS is
429    stored into that result. */
430 void
431 moments1_calculate (const struct moments1 *m,
432                     double *weight,
433                     double *mean, double *variance,
434                     double *skewness, double *kurtosis) 
435 {
436   assert (m != NULL);
437
438   if (mean != NULL)
439     *mean = SYSMIS;
440   if (variance != NULL)
441     *variance = SYSMIS;
442   if (skewness != NULL)
443     *skewness = SYSMIS;
444   if (kurtosis != NULL)
445     *kurtosis = SYSMIS;
446
447   if (weight != NULL)
448     *weight = m->w;
449
450   if (m->w > 0.) 
451     {
452       if (mean != NULL)
453         *mean = m->d1;
454
455       calc_moments (m->max_moment,
456                     m->w, 0., m->d2, m->d3, m->d4,
457                     variance, skewness, kurtosis);
458     }
459 }
460
461 /* Destroy one-pass moments M. */
462 void
463 moments1_destroy (struct moments1 *m) 
464 {
465   free (m);
466 }
467 \f
468 /* Returns the standard error of the skewness for the given total
469    weight W.
470
471    From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
472    section "DESCRIPTIVES". */
473 double
474 calc_seskew (double W)
475 {
476   return sqrt ((6. * W * (W - 1.)) / ((W - 2.) * (W + 1.) * (W + 3.)));
477 }
478
479 /* Returns the standard error of the kurtosis for the given total
480    weight W.
481
482    From _SPSS Statistical Algorithms, 2nd ed., 0-918469-89-9,
483    section "DESCRIPTIVES", except that the sqrt symbol is omitted
484    there. */
485 double
486 calc_sekurt (double W)
487 {
488   return sqrt ((4. * (pow2 (W) - 1.) * pow2 (calc_seskew (W)))
489                / ((W - 3.) * (W + 5.)));
490 }
491 \f
492 #include <stdio.h>
493 #include "command.h"
494 #include "lexer.h"
495
496 static int
497 read_values (double **values, double **weights, size_t *cnt) 
498 {
499   size_t cap = 0;
500
501   *values = NULL;
502   *weights = NULL;
503   *cnt = 0;
504   while (token == T_NUM) 
505     {
506       double value = tokval;
507       double weight = 1.;
508       lex_get ();
509       if (lex_match ('*'))
510         {
511           if (token != T_NUM) 
512             {
513               lex_error (_("expecting weight value"));
514               return 0;
515             }
516           weight = tokval;
517           lex_get ();
518         }
519
520       if (*cnt >= cap) 
521         {
522           cap = 2 * (cap + 8);
523           *values = xrealloc (*values, sizeof **values * cap);
524           *weights = xrealloc (*weights, sizeof **weights * cap);
525         }
526
527       (*values)[*cnt] = value;
528       (*weights)[*cnt] = weight;
529       (*cnt)++;
530     }
531
532   return 1;
533 }
534
535 int
536 cmd_debug_moments (void) 
537 {
538   int retval = CMD_FAILURE;
539   double *values = NULL;
540   double *weights = NULL;
541   double weight, M[4];
542   int two_pass = 1;
543   size_t cnt;
544   size_t i;
545
546   if (lex_match_id ("ONEPASS"))
547     two_pass = 0;
548   if (token != '/') 
549     {
550       lex_force_match ('/');
551       goto done;
552     }
553   fprintf (stderr, "%s => ", lex_rest_of_line (NULL));
554   lex_get ();
555
556   if (two_pass) 
557     {
558       struct moments *m = NULL;
559   
560       m = moments_create (MOMENT_KURTOSIS);
561       if (!read_values (&values, &weights, &cnt)) 
562         {
563           moments_destroy (m);
564           goto done; 
565         }
566       for (i = 0; i < cnt; i++)
567         moments_pass_one (m, values[i], weights[i]); 
568       for (i = 0; i < cnt; i++)
569         moments_pass_two (m, values[i], weights[i]);
570       moments_calculate (m, &weight, &M[0], &M[1], &M[2], &M[3]);
571       moments_destroy (m);
572     }
573   else 
574     {
575       struct moments1 *m = NULL;
576   
577       m = moments1_create (MOMENT_KURTOSIS);
578       if (!read_values (&values, &weights, &cnt)) 
579         {
580           moments1_destroy (m);
581           goto done; 
582         }
583       for (i = 0; i < cnt; i++)
584         moments1_add (m, values[i], weights[i]);
585       moments1_calculate (m, &weight, &M[0], &M[1], &M[2], &M[3]);
586       moments1_destroy (m);
587     }
588   
589   fprintf (stderr, "W=%.3f", weight);
590   for (i = 0; i < 4; i++) 
591     {
592       fprintf (stderr, " M%d=", i + 1);
593       if (M[i] == SYSMIS)
594         fprintf (stderr, "sysmis");
595       else if (fabs (M[i]) <= 0.0005)
596         fprintf (stderr, "0.000");
597       else
598         fprintf (stderr, "%.3f", M[i]);
599     }
600   fprintf (stderr, "\n");
601
602   retval = lex_end_of_command ();
603   
604  done:
605   free (values);
606   free (weights);
607   return retval;
608 }