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