src/math/levene.c

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