1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2010 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/>. */
23 #include <libpspp/misc.h>
25 #include <data/case.h>
26 #include <data/casereader.h>
27 #include <data/variable.h>
30 #include <libpspp/hmap.h>
34 struct hmap_node node;
45 find_group (struct hmap *map, const union value *target, int width)
48 unsigned int hash = value_hash (target, width, 0);
49 HMAP_FOR_EACH_WITH_HASH (l, struct lev, node, hash, map)
51 if (value_equal (&l->group, target, width))
61 levene (struct casereader *rx, const struct variable *gvar,
62 const struct variable *var, const struct variable *wv,
63 enum mv_class exclude)
65 double numerator = 0.0;
66 double denominator = 0.0;
68 double z_grand_mean = 0.0;
71 struct hmap map = HMAP_INITIALIZER (map);
74 struct casereader *r = casereader_clone (rx);
76 for (; (c = casereader_read (r)) != NULL; case_unref (c))
79 const union value *target = case_data (c, gvar);
80 int width = var_get_width (gvar);
81 unsigned int hash = value_hash (target, width, 0);
82 const double x = case_data (c, var)->f;
83 const double weight = wv ? case_data (c, wv)->f : 1.0;
85 if (var_is_value_missing (var, case_data (c, var), exclude))
88 l = find_group (&map, target, width);
92 l = xzalloc (sizeof *l);
93 value_clone (&l->group, target, width);
94 hmap_insert (&map, &l->node, hash);
98 l->t_bar += x * weight;
101 casereader_destroy (r);
105 HMAP_FOR_EACH (l, struct lev, node, &map)
111 n_groups = hmap_count (&map);
113 r = casereader_clone (rx);
114 for (; (c = casereader_read (r)) != NULL; case_unref (c))
116 struct lev *l = NULL;
117 const union value *target = case_data (c, gvar);
118 int width = var_get_width (gvar);
119 const double x = case_data (c, var)->f;
120 const double weight = wv ? case_data (c, wv)->f : 1.0;
122 if (var_is_value_missing (var, case_data (c, var), exclude))
125 l = find_group (&map, target, width);
128 l->z_mean += fabs (x - l->t_bar) * weight;
129 z_grand_mean += fabs (x - l->t_bar) * weight;
131 casereader_destroy (r);
135 HMAP_FOR_EACH (l, struct lev, node, &map)
141 z_grand_mean /= grand_n;
143 r = casereader_clone (rx);
144 for (; (c = casereader_read (r)) != NULL; case_unref (c))
148 const union value *target = case_data (c, gvar);
149 int width = var_get_width (gvar);
150 const double x = case_data (c, var)->f;
151 const double weight = wv ? case_data (c, wv)->f : 1.0;
153 if (var_is_value_missing (var, case_data (c, var), exclude))
156 l = find_group (&map, target, width);
159 z = fabs (x - l->t_bar);
160 denominator += pow2 (z - l->z_mean) * weight;
162 casereader_destroy (r);
164 denominator *= n_groups - 1;
167 double grand_n = 0.0;
170 HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &map)
172 numerator += l->n * pow2 (l->z_mean - z_grand_mean);
175 /* We don't need these anymore */
178 numerator *= grand_n - n_groups ;
182 return numerator/ denominator;