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