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       const union value *target = case_data (c, gvar);
 117       int width = var_get_width (gvar);
 118       const double x = case_data (c, var)->f;
 119       const double weight = wv ? case_data (c, wv)->f : 1.0;
 120
 121       if (var_is_value_missing (var, case_data (c, var), exclude))
 122         continue;
 123
 124       struct lev *l = find_group (&map, target, width);
 125
 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       const union value *target = case_data (c, gvar);
 148       int width = var_get_width (gvar);
 149       const double x = case_data (c, var)->f;
 150       const double weight = wv ? case_data (c, wv)->f : 1.0;
 151
 152       if (var_is_value_missing (var, case_data (c, var), exclude))
 153         continue;
 154
 155       struct lev *l = find_group (&map, target, width);
 156       assert (l);
 157       z = fabs (x - l->t_bar);
 158
 159       denominator += pow2 (z - l->z_mean) * weight;
 160     }
 161   casereader_destroy (r);
 162
 163   denominator *= n_groups - 1;
 164
 165   {
 166     double grand_n = 0.0;
 167     struct hmap_node *next;
 168     struct lev *l;
 169     HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &map)
 170       {
 171         numerator += l->n * pow2 (l->z_mean - z_grand_mean);
 172         grand_n += l->n;
 173
 174         /* We don't need these anymore */
 175         free (l);
 176       }
 177     numerator *= grand_n - n_groups ;
 178     hmap_destroy (&map);
 179   }
 180
 181   return numerator/ denominator;
 182 }