pintos-os.org Git - pspp/blob - 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 "levene.h"
  20 #include <math.h>
  21
  22 #include "libpspp/misc.h"
  23 #include "libpspp/hmap.h"
  24 #include "data/value.h"
  25 #include "data/val-type.h"
  26
  27 #include <gl/xalloc.h>
  28 #include <assert.h>
  29
  30 struct lev
  31 {
  32   struct hmap_node node;
  33   union value group;
  34
  35   double t_bar;
  36   double z_mean;
  37   double n;
  38 };
  39
  40 typedef unsigned int hash_func (const struct levene *, const union value *v);
  41 typedef bool cmp_func (const struct levene *, const union value *v0, const union value *v1);
  42
  43 struct levene
  44 {
  45   /* Width of the categorical variable */
  46   int gvw ;
  47
  48   /* The value deviding the the groups. Valid only for dichotomous categorical variable.*/
  49   const union value *cutpoint;
  50
  51
  52   /* A hashtable of struct lev objects indexed by union value */
  53   struct hmap hmap;
  54
  55   hash_func *hash;
  56   cmp_func *cmp;
  57
  58
  59   /* A state variable indicating how many passes have been done */
  60   int pass;
  61
  62   double grand_n;
  63   double z_grand_mean;
  64
  65   double denominator;
  66 };
  67
  68
  69 static unsigned int
  70 unique_hash (const struct levene *nl, const union value *val)
  71 {
  72   return value_hash (val, nl->gvw, 0);
  73 }
  74
  75 static bool
  76 unique_cmp (const struct levene *nl, const union value *val0, const union value *val1)
  77 {
  78   return value_equal (val0, val1, nl->gvw);
  79 }
  80
  81 static unsigned int
  82 cutpoint_hash (const struct levene *nl, const union value *val)
  83 {
  84   int x = value_compare_3way (val, nl->cutpoint, nl->gvw);
  85
  86   return (x < 0);
  87 }
  88
  89 static bool
  90 cutpoint_cmp (const struct levene *nl, const union value *val0, const union value *val1)
  91 {
  92   int x = value_compare_3way (val0, nl->cutpoint, nl->gvw);
  93
  94   int y = value_compare_3way (val1, nl->cutpoint, nl->gvw);
  95
  96   if ( x == 0) x = 1;
  97   if ( y == 0) y = 1;
  98
  99   return ( x == y);
 100 }
 101
 102
 103
 104 static struct lev *
 105 find_group (const struct levene *nl, const union value *target)
 106 {
 107   struct lev *l = NULL;
 108
 109   HMAP_FOR_EACH_WITH_HASH (l, struct lev, node, nl->hash (nl, target), &nl->hmap)
 110     {
 111       if (nl->cmp (nl, &l->group, target))
 112         break;
 113       l = NULL;
 114     }
 115   return l;
 116 }
 117
 118
 119 struct levene *
 120 levene_create (int indep_width, const union value *cutpoint)
 121 {
 122   struct levene *nl = xzalloc (sizeof *nl);
 123
 124   hmap_init (&nl->hmap);
 125
 126   nl->gvw = indep_width;
 127   nl->cutpoint = cutpoint;
 128
 129   nl->hash  = cutpoint ? cutpoint_hash : unique_hash;
 130   nl->cmp   = cutpoint ? cutpoint_cmp : unique_cmp;
 131
 132   return nl;
 133 }
 134
 135
 136 /* Data accumulation. First pass */
 137 void
 138 levene_pass_one (struct levene *nl, double value, double weight, const union value *gv)
 139 {
 140   struct lev *lev = find_group (nl, gv);
 141
 142   if ( nl->pass == 0 )
 143     {
 144       nl->pass = 1;
 145     }
 146   assert (nl->pass == 1);
 147
 148   if ( NULL == lev)
 149     {
 150       struct lev *l = xzalloc (sizeof *l);
 151       value_clone (&l->group, gv, nl->gvw);
 152       hmap_insert (&nl->hmap, &l->node, nl->hash (nl, &l->group));
 153       lev = l;
 154     }
 155
 156   lev->n += weight;
 157   lev->t_bar += value * weight;
 158
 159   nl->grand_n += weight;
 160 }
 161
 162 /* Data accumulation. Second pass */
 163 void
 164 levene_pass_two (struct levene *nl, double value, double weight, const union value *gv)
 165 {
 166   struct lev *lev = NULL;
 167
 168   if ( nl->pass == 1 )
 169     {
 170       struct lev *next;
 171       struct lev *l;
 172
 173       nl->pass = 2;
 174
 175       HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
 176       {
 177         l->t_bar /= l->n;
 178       }
 179     }
 180   assert (nl->pass == 2);
 181
 182   lev = find_group (nl, gv);
 183
 184   lev->z_mean += fabs (value - lev->t_bar) * weight;
 185   nl->z_grand_mean += fabs (value - lev->t_bar) * weight;
 186 }
 187
 188 /* Data accumulation. Third pass */
 189 void
 190 levene_pass_three (struct levene *nl, double value, double weight, const union value *gv)
 191 {
 192   double z;
 193   struct lev *lev = NULL;
 194
 195   if ( nl->pass == 2 )
 196     {
 197       struct lev *next;
 198       struct lev *l;
 199
 200       nl->pass = 3;
 201
 202       HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
 203       {
 204         l->z_mean /= l->n;
 205       }
 206
 207       nl->z_grand_mean /= nl->grand_n;
 208   }
 209
 210   assert (nl->pass == 3);
 211   lev = find_group (nl, gv);
 212
 213   z = fabs (value - lev->t_bar);
 214   nl->denominator += pow2 (z - lev->z_mean) * weight;
 215 }
 216
 217
 218 /* Return the value of the levene statistic */
 219 double
 220 levene_calculate (struct levene *nl)
 221 {
 222   struct lev *next;
 223   struct lev *l;
 224
 225   double numerator = 0.0;
 226   double nn = 0.0;
 227
 228   /* The Levene calculation requires three passes.
 229      Normally this should have been done prior to calling this function.
 230      However, in abnormal circumstances (eg. the dataset is empty) there
 231      will have been no passes.
 232    */
 233   assert (nl->pass == 0 || nl->pass == 3);
 234
 235   if ( nl->pass == 0 )
 236     return SYSMIS;
 237
 238   nl->denominator *= hmap_count (&nl->hmap) - 1;
 239
 240   HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
 241     {
 242       numerator += l->n * pow2 (l->z_mean - nl->z_grand_mean);
 243       nn += l->n;
 244     }
 245
 246   numerator *= nn - hmap_count (&nl->hmap);
 247
 248   return numerator / nl->denominator;
 249 }
 250
 251 void
 252 levene_destroy (struct levene *nl)
 253 {
 254   struct lev *next;
 255   struct lev *l;
 256
 257   HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
 258     {
 259       value_destroy (&l->group, nl->gvw);
 260       free (l);
 261     }
 262
 263   hmap_destroy (&nl->hmap);
 264   free (nl);
 265 }