gui: Fix behavior of print preview in psppire-output-window.
[pspp] / src / math / levene.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 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 "levene.h"
20 #include <math.h>
21
22 #include "libpspp/misc.h"
23 #include "libpspp/hmap.h"
24 #include "data/value.h"
25
26 #include <gl/xalloc.h>
27 #include <assert.h>
28
29 struct lev
30 {
31   struct hmap_node node;
32   union value group;
33
34   double t_bar;
35   double z_mean;
36   double n;
37 };
38
39 typedef unsigned int hash_func (const struct levene *, const union value *v);
40 typedef bool cmp_func (const struct levene *, const union value *v0, const union value *v1);
41
42 struct levene
43 {
44   /* Width of the categorical variable */
45   int gvw ;
46
47   /* The value deviding the the groups. Valid only for dichotomous categorical variable.*/
48   const union value *cutpoint;
49
50
51   /* A hashtable of struct lev objects indexed by union value */
52   struct hmap hmap;
53
54   hash_func *hash;
55   cmp_func *cmp;
56
57
58   /* A state variable indicating how many passes have been done */
59   int pass;
60
61   double grand_n;
62   double z_grand_mean;
63
64   double denominator;
65 };
66
67
68 static unsigned int
69 unique_hash (const struct levene *nl, const union value *val)
70 {
71   return value_hash (val, nl->gvw, 0);
72 }
73
74 static bool
75 unique_cmp (const struct levene *nl, const union value *val0, const union value *val1)
76 {
77   return value_equal (val0, val1, nl->gvw);
78 }
79
80 static unsigned int
81 cutpoint_hash (const struct levene *nl, const union value *val)
82 {
83   int x = value_compare_3way (val, nl->cutpoint, nl->gvw);
84
85   return (x < 0);
86 }
87
88 static bool
89 cutpoint_cmp (const struct levene *nl, const union value *val0, const union value *val1)
90 {
91   int x = value_compare_3way (val0, nl->cutpoint, nl->gvw);
92
93   int y = value_compare_3way (val1, nl->cutpoint, nl->gvw);
94
95   if ( x == 0) x = 1;
96   if ( y == 0) y = 1;
97
98   return ( x == y);
99 }
100
101
102
103 static struct lev *
104 find_group (const struct levene *nl, const union value *target)
105 {
106   struct lev *l = NULL;
107
108   HMAP_FOR_EACH_WITH_HASH (l, struct lev, node, nl->hash (nl, target), &nl->hmap)
109     {
110       if (nl->cmp (nl, &l->group, target))
111         break;
112       l = NULL;
113     }
114   return l;
115 }
116
117
118 struct levene *
119 levene_create (int indep_width, const union value *cutpoint)
120 {
121   struct levene *nl = xzalloc (sizeof *nl);
122
123   hmap_init (&nl->hmap);
124
125   nl->gvw = indep_width;
126   nl->cutpoint = cutpoint;
127
128   nl->hash  = cutpoint ? cutpoint_hash : unique_hash;
129   nl->cmp   = cutpoint ? cutpoint_cmp : unique_cmp;
130
131   return nl;
132 }
133
134
135 /* Data accumulation. First pass */
136 void 
137 levene_pass_one (struct levene *nl, double value, double weight, const union value *gv)
138 {
139   struct lev *lev = find_group (nl, gv);
140
141   if ( nl->pass == 0 ) 
142     {
143       nl->pass = 1;
144     }
145   assert (nl->pass == 1);
146
147   if ( NULL == lev)
148     {
149       struct lev *l = xzalloc (sizeof *l);
150       value_clone (&l->group, gv, nl->gvw);
151       hmap_insert (&nl->hmap, &l->node, nl->hash (nl, &l->group));
152       lev = l;
153     }
154
155   lev->n += weight;
156   lev->t_bar += value * weight;
157
158   nl->grand_n += weight;
159 }
160
161 /* Data accumulation. Second pass */
162 void 
163 levene_pass_two (struct levene *nl, double value, double weight, const union value *gv)
164 {
165   struct lev *lev = NULL;
166
167   if ( nl->pass == 1 )
168     {
169       struct lev *next;
170       struct lev *l;
171
172       nl->pass = 2;
173
174       HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
175       {
176         l->t_bar /= l->n;
177       }
178     }
179   assert (nl->pass == 2);
180
181   lev = find_group (nl, gv);
182
183   lev->z_mean += fabs (value - lev->t_bar) * weight;
184   nl->z_grand_mean += fabs (value - lev->t_bar) * weight;
185 }
186
187 /* Data accumulation. Third pass */
188 void 
189 levene_pass_three (struct levene *nl, double value, double weight, const union value *gv)
190 {
191   double z;
192   struct lev *lev = NULL;
193
194   if ( nl->pass == 2 )
195     {
196       struct lev *next;
197       struct lev *l;
198
199       nl->pass = 3;
200
201       HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
202       {
203         l->z_mean /= l->n;
204       }
205
206       nl->z_grand_mean /= nl->grand_n;
207   }
208
209   assert (nl->pass == 3);
210   lev = find_group (nl, gv);
211
212   z = fabs (value - lev->t_bar);
213   nl->denominator += pow2 (z - lev->z_mean) * weight;
214 }
215
216
217 /* Return the value of the levene statistic */
218 double
219 levene_calculate (struct levene *nl)
220 {
221   struct lev *next;
222   struct lev *l;
223
224   double numerator = 0.0;
225   double nn = 0.0;
226   if ( nl->pass == 3 )
227     {
228       nl->pass = 4;
229     }
230
231   assert (nl->pass == 4);
232
233   nl->denominator *= hmap_count (&nl->hmap) - 1;
234
235   HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
236     {
237       numerator += l->n * pow2 (l->z_mean - nl->z_grand_mean);
238       nn += l->n;
239     }
240
241   numerator *= nn - hmap_count (&nl->hmap);
242     
243   return numerator / nl->denominator;
244 }
245
246 void
247 levene_destroy (struct levene *nl)
248 {
249   struct lev *next;
250   struct lev *l;
251
252   HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
253     {
254       value_destroy (&l->group, nl->gvw);
255       free (l);
256     }
257
258   hmap_destroy (&nl->hmap);
259   free (nl);
260 }