-/* This file is part of GNU PSPP
- Computes Levene test statistic.
+/* PSPP - a program for statistical analysis.
+ Copyright (C) 2010, 2011 Free Software Foundation, Inc.
- Copyright (C) 2004 Free Software Foundation, Inc.
- Written by John Darrington <john@darrington.wattle.id.au>
+ This program is free software: you can redistribute it and/or modify
+ it under the terms of the GNU General Public License as published by
+ the Free Software Foundation, either version 3 of the License, or
+ (at your option) any later version.
- This program is free software; you can redistribute it and/or
- modify it under the terms of the GNU General Public License as
- published by the Free Software Foundation; either version 2 of the
- License, or (at your option) any later version.
-
- This program is distributed in the hope that it will be useful, but
- WITHOUT ANY WARRANTY; without even the implied warranty of
- MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
- General Public License for more details.
+ This program is distributed in the hope that it will be useful,
+ but WITHOUT ANY WARRANTY; without even the implied warranty of
+ MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
+ GNU General Public License for more details.
You should have received a copy of the GNU General Public License
- along with this program; if not, write to the Free Software
- Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA
- 02110-1301, USA. */
+ along with this program. If not, see <http://www.gnu.org/licenses/>. */
#include <config.h>
-#include "levene.h"
-#include <libpspp/message.h>
-#include <data/case.h>
-#include <data/casefile.h>
-#include <data/dictionary.h>
-#include "group-proc.h"
-#include <libpspp/hash.h>
-#include <libpspp/str.h>
-#include <data/variable.h>
-#include <procedure.h>
-#include <libpspp/alloc.h>
-#include <libpspp/misc.h>
-#include "group.h"
+#include "levene.h"
#include <math.h>
-#include <stdlib.h>
-
-/* This module calculates the Levene statistic for variables.
+#include "libpspp/misc.h"
+#include "libpspp/hmap.h"
+#include "data/value.h"
- Just for reference, the Levene Statistic is a defines as follows:
+#include <gl/xalloc.h>
+#include <assert.h>
- W = \frac{ (n-k)\sum_{i=1}^k n_i(Z_{iL} - Z_{LL})^2}
- { (k-1)\sum_{i=1}^k \sum_{j=1}^{n_i} (Z_{ij} - Z_{iL})^2}
+struct lev
+{
+ struct hmap_node node;
+ union value group;
- where:
- k is the number of groups
- n is the total number of samples
- n_i is the number of samples in the ith group
- Z_{ij} is | Y_{ij} - Y_{iL} | where Y_{iL} is the mean of the ith group
- Z_{iL} is the mean of Z_{ij} over the ith group
- Z_{LL} is the grand mean of Z_{ij}
+ double t_bar;
+ double z_mean;
+ double n;
+};
- Imagine calculating that with pencil and paper!
+typedef unsigned int hash_func (const struct levene *, const union value *v);
+typedef bool cmp_func (const struct levene *, const union value *v0, const union value *v1);
- */
+struct levene
+{
+ /* Width of the categorical variable */
+ int gvw ;
+ /* The value deviding the the groups. Valid only for dichotomous categorical variable.*/
+ const union value *cutpoint;
-struct levene_info
-{
- /* Per group statistics */
- struct t_test_proc **group_stats;
+ /* A hashtable of struct lev objects indexed by union value */
+ struct hmap hmap;
- /* The independent variable */
- struct variable *v_indep;
+ hash_func *hash;
+ cmp_func *cmp;
- /* Number of dependent variables */
- size_t n_dep;
- /* The dependent variables */
- struct variable **v_dep;
+ /* A state variable indicating how many passes have been done */
+ int pass;
- /* How to treat missing values */
- enum lev_missing missing;
+ double grand_n;
+ double z_grand_mean;
- /* Function to test for missing values */
- is_missing_func *is_missing;
+ double denominator;
};
-/* First pass */
-static void levene_precalc (const struct levene_info *l);
-static int levene_calc (const struct ccase *, void *);
-static void levene_postcalc (void *);
+static unsigned int
+unique_hash (const struct levene *nl, const union value *val)
+{
+ return value_hash (val, nl->gvw, 0);
+}
-/* Second pass */
-static void levene2_precalc (void *);
-static int levene2_calc (const struct ccase *, void *);
-static void levene2_postcalc (void *);
+static bool
+unique_cmp (const struct levene *nl, const union value *val0, const union value *val1)
+{
+ return value_equal (val0, val1, nl->gvw);
+}
+static unsigned int
+cutpoint_hash (const struct levene *nl, const union value *val)
+{
+ int x = value_compare_3way (val, nl->cutpoint, nl->gvw);
+
+ return (x < 0);
+}
-void
-levene(const struct casefile *cf,
- struct variable *v_indep, size_t n_dep, struct variable **v_dep,
- enum lev_missing missing, is_missing_func value_is_missing)
+static bool
+cutpoint_cmp (const struct levene *nl, const union value *val0, const union value *val1)
{
- struct casereader *r;
- struct ccase c;
- struct levene_info l;
+ int x = value_compare_3way (val0, nl->cutpoint, nl->gvw);
- l.n_dep = n_dep;
- l.v_indep = v_indep;
- l.v_dep = v_dep;
- l.missing = missing;
- l.is_missing = value_is_missing;
+ int y = value_compare_3way (val1, nl->cutpoint, nl->gvw);
+ if ( x == 0) x = 1;
+ if ( y == 0) y = 1;
+
+ return ( x == y);
+}
- levene_precalc(&l);
- for(r = casefile_get_reader (cf);
- casereader_read (r, &c) ;
- case_destroy (&c))
- {
- levene_calc(&c,&l);
- }
- casereader_destroy (r);
- levene_postcalc(&l);
- levene2_precalc(&l);
- for(r = casefile_get_reader (cf);
- casereader_read (r, &c) ;
- case_destroy (&c))
+static struct lev *
+find_group (const struct levene *nl, const union value *target)
+{
+ struct lev *l = NULL;
+
+ HMAP_FOR_EACH_WITH_HASH (l, struct lev, node, nl->hash (nl, target), &nl->hmap)
{
- levene2_calc(&c,&l);
+ if (nl->cmp (nl, &l->group, target))
+ break;
+ l = NULL;
}
- casereader_destroy (r);
- levene2_postcalc(&l);
-
+ return l;
}
-/* Internal variables used in calculating the Levene statistic */
-/* Per variable statistics */
-struct lz_stats
+struct levene *
+levene_create (int indep_width, const union value *cutpoint)
{
- /* Total of all lz */
- double grand_total;
+ struct levene *nl = xzalloc (sizeof *nl);
- /* Mean of all lz */
- double grand_mean;
+ hmap_init (&nl->hmap);
- /* The total number of cases */
- double total_n ;
+ nl->gvw = indep_width;
+ nl->cutpoint = cutpoint;
- /* Number of groups */
- int n_groups;
-};
+ nl->hash = cutpoint ? cutpoint_hash : unique_hash;
+ nl->cmp = cutpoint ? cutpoint_cmp : unique_cmp;
-/* An array of lz_stats for each variable */
-static struct lz_stats *lz;
+ return nl;
+}
-static void
-levene_precalc (const struct levene_info *l)
+/* Data accumulation. First pass */
+void
+levene_pass_one (struct levene *nl, double value, double weight, const union value *gv)
{
- size_t i;
-
- lz = xnmalloc (l->n_dep, sizeof *lz);
+ struct lev *lev = find_group (nl, gv);
- for(i = 0; i < l->n_dep ; ++i )
+ if ( nl->pass == 0 )
{
- struct variable *var = l->v_dep[i];
- struct group_proc *gp = group_proc_get (var);
- struct group_statistics *gs;
- struct hsh_iterator hi;
-
- lz[i].grand_total = 0;
- lz[i].total_n = 0;
- lz[i].n_groups = gp->n_groups ;
-
-
- for ( gs = hsh_first(gp->group_hash, &hi);
- gs != 0;
- gs = hsh_next(gp->group_hash, &hi))
- {
- gs->lz_total = 0;
- }
-
+ nl->pass = 1;
}
+ assert (nl->pass == 1);
-}
-
-static int
-levene_calc (const struct ccase *c, void *_l)
-{
- size_t i;
- int warn = 0;
- struct levene_info *l = (struct levene_info *) _l;
- const union value *gv = case_data (c, l->v_indep->fv);
- struct group_statistics key;
- double weight = dict_get_case_weight(default_dict,c,&warn);
-
- /* Skip the entire case if /MISSING=LISTWISE is set */
- if ( l->missing == LEV_LISTWISE )
+ if ( NULL == lev)
{
- for (i = 0; i < l->n_dep; ++i)
- {
- struct variable *v = l->v_dep[i];
- const union value *val = case_data (c, v->fv);
-
- if (l->is_missing (&v->miss, val) )
- {
- return 0;
- }
- }
+ struct lev *l = xzalloc (sizeof *l);
+ value_clone (&l->group, gv, nl->gvw);
+ hmap_insert (&nl->hmap, &l->node, nl->hash (nl, &l->group));
+ lev = l;
}
-
- key.id = *gv;
-
- for (i = 0; i < l->n_dep; ++i)
- {
- struct variable *var = l->v_dep[i];
- struct group_proc *gp = group_proc_get (var);
- double levene_z;
- const union value *v = case_data (c, var->fv);
- struct group_statistics *gs;
-
- gs = hsh_find(gp->group_hash,(void *) &key );
-
- if ( 0 == gs )
- continue ;
-
- if ( ! l->is_missing(&var->miss, v))
- {
- levene_z= fabs(v->f - gs->mean);
- lz[i].grand_total += levene_z * weight;
- lz[i].total_n += weight;
-
- gs->lz_total += levene_z * weight;
- }
+ lev->n += weight;
+ lev->t_bar += value * weight;
- }
- return 0;
+ nl->grand_n += weight;
}
-
-static void
-levene_postcalc (void *_l)
+/* Data accumulation. Second pass */
+void
+levene_pass_two (struct levene *nl, double value, double weight, const union value *gv)
{
- size_t v;
+ struct lev *lev = NULL;
- struct levene_info *l = (struct levene_info *) _l;
-
- for (v = 0; v < l->n_dep; ++v)
+ if ( nl->pass == 1 )
{
- /* This is Z_LL */
- lz[v].grand_mean = lz[v].grand_total / lz[v].total_n ;
- }
-
-
-}
+ struct lev *next;
+ struct lev *l;
+ nl->pass = 2;
-/* The denominator for the expression for the Levene */
-static double *lz_denominator;
+ HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
+ {
+ l->t_bar /= l->n;
+ }
+ }
+ assert (nl->pass == 2);
-static void
-levene2_precalc (void *_l)
-{
- size_t v;
+ lev = find_group (nl, gv);
- struct levene_info *l = (struct levene_info *) _l;
+ lev->z_mean += fabs (value - lev->t_bar) * weight;
+ nl->z_grand_mean += fabs (value - lev->t_bar) * weight;
+}
- lz_denominator = xnmalloc (l->n_dep, sizeof *lz_denominator);
+/* Data accumulation. Third pass */
+void
+levene_pass_three (struct levene *nl, double value, double weight, const union value *gv)
+{
+ double z;
+ struct lev *lev = NULL;
- /* This stuff could go in the first post calc . . . */
- for (v = 0; v < l->n_dep; ++v)
+ if ( nl->pass == 2 )
{
- struct hsh_iterator hi;
- struct group_statistics *g;
+ struct lev *next;
+ struct lev *l;
- struct variable *var = l->v_dep[v] ;
- struct hsh_table *hash = group_proc_get (var)->group_hash;
+ nl->pass = 3;
+ HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
+ {
+ l->z_mean /= l->n;
+ }
- for(g = (struct group_statistics *) hsh_first(hash,&hi);
- g != 0 ;
- g = (struct group_statistics *) hsh_next(hash,&hi) )
- {
- g->lz_mean = g->lz_total / g->n ;
- }
- lz_denominator[v] = 0;
+ nl->z_grand_mean /= nl->grand_n;
}
-}
-static int
-levene2_calc (const struct ccase *c, void *_l)
-{
- size_t i;
- int warn = 0;
+ assert (nl->pass == 3);
+ lev = find_group (nl, gv);
- struct levene_info *l = (struct levene_info *) _l;
+ z = fabs (value - lev->t_bar);
+ nl->denominator += pow2 (z - lev->z_mean) * weight;
+}
- double weight = dict_get_case_weight(default_dict,c,&warn);
- const union value *gv = case_data (c, l->v_indep->fv);
- struct group_statistics key;
+/* Return the value of the levene statistic */
+double
+levene_calculate (struct levene *nl)
+{
+ struct lev *next;
+ struct lev *l;
- /* Skip the entire case if /MISSING=LISTWISE is set */
- if ( l->missing == LEV_LISTWISE )
+ double numerator = 0.0;
+ double nn = 0.0;
+ if ( nl->pass == 3 )
{
- for (i = 0; i < l->n_dep; ++i)
- {
- struct variable *v = l->v_dep[i];
- const union value *val = case_data (c, v->fv);
-
- if (l->is_missing(&v->miss, val) )
- {
- return 0;
- }
- }
+ nl->pass = 4;
}
- key.id = *gv;
-
- for (i = 0; i < l->n_dep; ++i)
- {
- double levene_z;
- struct variable *var = l->v_dep[i] ;
- const union value *v = case_data (c, var->fv);
- struct group_statistics *gs;
-
- gs = hsh_find(group_proc_get (var)->group_hash,(void *) &key );
+ assert (nl->pass == 4);
- if ( 0 == gs )
- continue;
+ nl->denominator *= hmap_count (&nl->hmap) - 1;
- if ( ! l->is_missing (&var->miss, v) )
- {
- levene_z = fabs(v->f - gs->mean);
- lz_denominator[i] += weight * pow2(levene_z - gs->lz_mean);
- }
+ HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
+ {
+ numerator += l->n * pow2 (l->z_mean - nl->z_grand_mean);
+ nn += l->n;
}
- return 0;
+ numerator *= nn - hmap_count (&nl->hmap);
+
+ return numerator / nl->denominator;
}
-
-static void
-levene2_postcalc (void *_l)
+void
+levene_destroy (struct levene *nl)
{
- size_t v;
+ struct lev *next;
+ struct lev *l;
- struct levene_info *l = (struct levene_info *) _l;
-
- for (v = 0; v < l->n_dep; ++v)
+ HMAP_FOR_EACH_SAFE (l, next, struct lev, node, &nl->hmap)
{
- double lz_numerator = 0;
- struct hsh_iterator hi;
- struct group_statistics *g;
-
- struct variable *var = l->v_dep[v] ;
- struct group_proc *gp = group_proc_get (var);
- struct hsh_table *hash = gp->group_hash;
-
- for(g = (struct group_statistics *) hsh_first(hash,&hi);
- g != 0 ;
- g = (struct group_statistics *) hsh_next(hash,&hi) )
- {
- lz_numerator += g->n * pow2(g->lz_mean - lz[v].grand_mean );
- }
- lz_numerator *= ( gp->ugs.n - gp->n_groups );
-
- lz_denominator[v] *= (gp->n_groups - 1);
-
- gp->levene = lz_numerator / lz_denominator[v] ;
-
+ value_destroy (&l->group, nl->gvw);
+ free (l);
}
- /* Now clear up after ourselves */
- free(lz_denominator);
- free(lz);
+ hmap_destroy (&nl->hmap);
+ free (nl);
}
-
projects / pspp / blobdiff