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 "math/levene.h"
  20
  21 #include <math.h>
  22
  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"
  28
  29 struct lev
  30 {
  31   struct hmap_node node;
  32   union value group;
  33   int width;
  34
  35   double t_bar;
  36   double z_mean;
  37   double n;
  38 };
  39
  40
  41 static struct lev *
  42 find_group (struct hmap *map, const union value *target, int width)
  43 {
  44   struct lev *l = NULL;
  45   unsigned int hash = value_hash (target, width, 0);
  46   HMAP_FOR_EACH_WITH_HASH (l, struct lev, node, hash, map)
  47     {
  48       if (value_equal (&l->group, target, width))
  49         break;
  50       l = NULL;
  51     }
  52
  53   return l;
  54 }
  55
  56
  57 double
  58 levene (struct casereader *rx, const struct variable *gvar,
  59         const struct variable *var, const struct variable *wv,
  60         enum mv_class exclude)
  61 {
  62   double numerator = 0.0;
  63   double denominator = 0.0;
  64   int n_groups = 0;
  65   double z_grand_mean = 0.0;
  66   double grand_n = 0.0;
  67
  68   struct hmap map = HMAP_INITIALIZER (map);
  69
  70   struct ccase *c;
  71   struct casereader *r = casereader_clone (rx);
  72
  73   for (; (c = casereader_read (r)) != NULL; case_unref (c))
  74     {
  75       struct lev *l = NULL;
  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;
  81
  82       if (var_is_value_missing (var, case_data (c, var), exclude))
  83         continue;
  84
  85       l = find_group (&map, target, width);
  86
  87       if (l == NULL)
  88         {
  89           l = xzalloc (sizeof *l);
  90           value_clone (&l->group, target, width);
  91           hmap_insert (&map, &l->node, hash);
  92         }
  93
  94       l->n += weight;
  95       l->t_bar += x * weight;
  96       grand_n += weight;
  97     }
  98   casereader_destroy (r);
  99
 100   {
 101     struct lev *l;
 102     HMAP_FOR_EACH (l, struct lev, node, &map)
 103       {
 104         l->t_bar /= l->n;
 105       }
 106   }
 107
 108   n_groups = hmap_count (&map);
 109
 110   r = casereader_clone (rx);
 111   for (; (c = casereader_read (r)) != NULL; case_unref (c))
 112     {
 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;
 118
 119       if (var_is_value_missing (var, case_data (c, var), exclude))
 120         continue;
 121
 122       l = find_group (&map, target, width);
 123       assert (l);
 124
 125       l->z_mean += fabs (x - l->t_bar) * weight;
 126       z_grand_mean += fabs (x - l->t_bar) * weight;
 127     }
 128   casereader_destroy (r);
 129
 130   {
 131     struct lev *l;
 132     HMAP_FOR_EACH (l, struct lev, node, &map)
 133       {
 134         l->z_mean /= l->n;
 135       }
 136   }
 137
 138   z_grand_mean /= grand_n;
 139
 140   r = casereader_clone (rx);
 141   for (; (c = casereader_read (r)) != NULL; case_unref (c))
 142     {
 143       double z;
 144       struct lev *l;
 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;
 149
 150       if (var_is_value_missing (var, case_data (c, var), exclude))
 151         continue;
 152
 153       l = find_group (&map, target, width);
 154       assert (l);
 155
 156       z = fabs (x - l->t_bar);
 157       denominator += pow2 (z - l->z_mean) * weight;
 158     }
 159   casereader_destroy (r);
 160
 161   denominator *= n_groups - 1;
 162
 163   {
 164     double grand_n = 0.0;
 165     struct lev *next;
 166     struct lev *l;
 167     HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &map)
 168       {
 169         numerator += l->n * pow2 (l->z_mean - z_grand_mean);
 170         grand_n += l->n;
 171
 172         /* We don't need these anymore */
 173         free (l);
 174       }
 175     numerator *= grand_n - n_groups ;
 176     hmap_destroy (&map);
 177   }
 178
 179   return numerator/ denominator;
 180 }