1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2010, 2011 Free Software Foundation, Inc.
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.
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.
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/>. */
19 #include "math/levene.h"
23 #include "data/case.h"
24 #include "data/casereader.h"
25 #include "data/variable.h"
26 #include "libpspp/hmap.h"
27 #include "libpspp/misc.h"
31 struct hmap_node node;
42 find_group (struct hmap *map, const union value *target, int width)
45 unsigned int hash = value_hash (target, width, 0);
46 HMAP_FOR_EACH_WITH_HASH (l, struct lev, node, hash, map)
48 if (value_equal (&l->group, target, width))
58 levene (struct casereader *rx, const struct variable *gvar,
59 const struct variable *var, const struct variable *wv,
60 enum mv_class exclude)
62 double numerator = 0.0;
63 double denominator = 0.0;
65 double z_grand_mean = 0.0;
68 struct hmap map = HMAP_INITIALIZER (map);
71 struct casereader *r = casereader_clone (rx);
73 for (; (c = casereader_read (r)) != NULL; case_unref (c))
76 const union value *target = case_data (c, gvar);
77 int width = var_get_width (gvar);
78 unsigned int hash = value_hash (target, width, 0);
79 const double x = case_data (c, var)->f;
80 const double weight = wv ? case_data (c, wv)->f : 1.0;
82 if (var_is_value_missing (var, case_data (c, var), exclude))
85 l = find_group (&map, target, width);
89 l = xzalloc (sizeof *l);
90 value_clone (&l->group, target, width);
91 hmap_insert (&map, &l->node, hash);
95 l->t_bar += x * weight;
98 casereader_destroy (r);
102 HMAP_FOR_EACH (l, struct lev, node, &map)
108 n_groups = hmap_count (&map);
110 r = casereader_clone (rx);
111 for (; (c = casereader_read (r)) != NULL; case_unref (c))
113 struct lev *l = NULL;
114 const union value *target = case_data (c, gvar);
115 int width = var_get_width (gvar);
116 const double x = case_data (c, var)->f;
117 const double weight = wv ? case_data (c, wv)->f : 1.0;
119 if (var_is_value_missing (var, case_data (c, var), exclude))
122 l = find_group (&map, target, width);
125 l->z_mean += fabs (x - l->t_bar) * weight;
126 z_grand_mean += fabs (x - l->t_bar) * weight;
128 casereader_destroy (r);
132 HMAP_FOR_EACH (l, struct lev, node, &map)
138 z_grand_mean /= grand_n;
140 r = casereader_clone (rx);
141 for (; (c = casereader_read (r)) != NULL; case_unref (c))
145 const union value *target = case_data (c, gvar);
146 int width = var_get_width (gvar);
147 const double x = case_data (c, var)->f;
148 const double weight = wv ? case_data (c, wv)->f : 1.0;
150 if (var_is_value_missing (var, case_data (c, var), exclude))
153 l = find_group (&map, target, width);
156 z = fabs (x - l->t_bar);
157 denominator += pow2 (z - l->z_mean) * weight;
159 casereader_destroy (r);
161 denominator *= n_groups - 1;
164 double grand_n = 0.0;
167 HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &map)
169 numerator += l->n * pow2 (l->z_mean - z_grand_mean);
172 /* We don't need these anymore */
175 numerator *= grand_n - n_groups ;
179 return numerator/ denominator;