f5fc858dd4d2b130f0af647f14ad7483e92dc311
[pspp] / src / language / tests / moments-test.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 1997-9, 2000, 2010 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 #include <stdio.h>
19 #include "gettext.h"
20 #include <language/command.h>
21 #include <language/lexer/lexer.h>
22 #include <math/moments.h>
23 #include <math.h>
24 #include <stdlib.h>
25 #include "xalloc.h"
26 #include <libpspp/compiler.h>
27
28 #define _(msgid) gettext (msgid)
29
30 static bool
31 read_values (struct lexer *lexer, double **values, double **weights, size_t *cnt)
32 {
33   size_t cap = 0;
34
35   *values = NULL;
36   *weights = NULL;
37   *cnt = 0;
38   while (lex_is_number (lexer))
39     {
40       double value = lex_tokval (lexer);
41       double weight = 1.;
42       lex_get (lexer);
43       if (lex_match (lexer, T_ASTERISK))
44         {
45           if (!lex_is_number (lexer))
46             {
47               lex_error (lexer, _("expecting weight value"));
48               return false;
49             }
50           weight = lex_tokval (lexer);
51           lex_get (lexer);
52         }
53
54       if (*cnt >= cap)
55         {
56           cap = 2 * (cap + 8);
57           *values = xnrealloc (*values, cap, sizeof **values);
58           *weights = xnrealloc (*weights, cap, sizeof **weights);
59         }
60
61       (*values)[*cnt] = value;
62       (*weights)[*cnt] = weight;
63       (*cnt)++;
64     }
65
66   return true;
67 }
68
69 int
70 cmd_debug_moments (struct lexer *lexer, struct dataset *ds UNUSED)
71 {
72   int retval = CMD_FAILURE;
73   double *values = NULL;
74   double *weights = NULL;
75   double weight, M[4];
76   int two_pass = 1;
77   size_t cnt;
78   size_t i;
79
80   if (lex_match_id (lexer, "ONEPASS"))
81     two_pass = 0;
82   if (lex_token (lexer) != T_SLASH)
83     {
84       lex_force_match (lexer, T_SLASH);
85       goto done;
86     }
87   lex_get (lexer);
88
89   if (two_pass)
90     {
91       struct moments *m = NULL;
92
93       m = moments_create (MOMENT_KURTOSIS);
94       if (!read_values (lexer, &values, &weights, &cnt))
95         {
96           moments_destroy (m);
97           goto done;
98         }
99       for (i = 0; i < cnt; i++)
100         moments_pass_one (m, values[i], weights[i]);
101       for (i = 0; i < cnt; i++)
102         moments_pass_two (m, values[i], weights[i]);
103       moments_calculate (m, &weight, &M[0], &M[1], &M[2], &M[3]);
104       moments_destroy (m);
105     }
106   else
107     {
108       struct moments1 *m = NULL;
109
110       m = moments1_create (MOMENT_KURTOSIS);
111       if (!read_values (lexer, &values, &weights, &cnt))
112         {
113           moments1_destroy (m);
114           goto done;
115         }
116       for (i = 0; i < cnt; i++)
117         moments1_add (m, values[i], weights[i]);
118       moments1_calculate (m, &weight, &M[0], &M[1], &M[2], &M[3]);
119       moments1_destroy (m);
120     }
121
122   fprintf (stderr, "W=%.3f", weight);
123   for (i = 0; i < 4; i++)
124     {
125       fprintf (stderr, " M%zu=", i + 1);
126       if (M[i] == SYSMIS)
127         fprintf (stderr, "sysmis");
128       else if (fabs (M[i]) <= 0.0005)
129         fprintf (stderr, "0.000");
130       else
131         fprintf (stderr, "%.3f", M[i]);
132     }
133   fprintf (stderr, "\n");
134
135   retval = lex_end_of_command (lexer);
136
137  done:
138   free (values);
139   free (weights);
140   return retval;
141 }