src/factor_stats.c

   1 /* PSPP - A program for statistical analysis . -*-c-*-
   2
   3 Copyright (C) 2004 Free Software Foundation, Inc.
   4 Author: John Darrington 2004
   5
   6 This program is free software; you can redistribute it and/or
   7 modify it under the terms of the GNU General Public License as
   8 published by the Free Software Foundation; either version 2 of the
   9 License, or (at your option) any later version.
  10
  11 This program is distributed in the hope that it will be useful, but
  12 WITHOUT ANY WARRANTY; without even the implied warranty of
  13 MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
  14 General Public License for more details.
  15
  16 You should have received a copy of the GNU General Public License
  17 along with this program; if not, write to the Free Software
  18 Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA
  19 02111-1307, USA. */
  20
  21 #include "factor_stats.h"
  22 #include "config.h"
  23 #include "val.h"
  24 #include "hash.h"
  25 #include "algorithm.h"
  26 #include "alloc.h"
  27
  28 #include <stdlib.h>
  29 #include <math.h>
  30 #include <float.h>
  31 #include <assert.h>
  32
  33
  34
  35 void
  36 metrics_precalc(struct metrics *fs)
  37 {
  38   fs->n = 0;
  39   fs->ssq = 0;
  40   fs->sum = 0;
  41   fs->min = DBL_MAX;
  42   fs->max = -DBL_MAX;
  43
  44   fs->ordered_data = hsh_create(20,
  45                                 (hsh_compare_func *) compare_values,
  46                                 (hsh_hash_func *) hash_value,
  47                                 0,
  48                                 (void *) 0);
  49 }
  50
  51 void
  52 metrics_calc(struct metrics *fs, const union value *val, double weight)
  53 {
  54
  55
  56   struct weighted_value **wv;
  57   const double x = val->f;
  58
  59   fs->n    += weight;
  60   fs->ssq  += x * x * weight;
  61   fs->sum  += x * weight;
  62
  63   if ( x < fs->min) fs->min = x;
  64   if ( x > fs->max) fs->max = x;
  65
  66
  67   wv = (struct weighted_value **) hsh_probe (fs->ordered_data,(void *) val );
  68
  69   if ( *wv  )
  70     {
  71       /* If this value has already been seen, then simply
  72          increase its weight */
  73
  74       assert( (*wv)->v.f == val->f );
  75       (*wv)->w += weight;
  76     }
  77   else
  78     {
  79       *wv = xmalloc( sizeof (struct weighted_value) );
  80       (*wv)->v = *val;
  81       (*wv)->w = weight;
  82       hsh_insert(fs->ordered_data,(void *) *wv);
  83     }
  84
  85 }
  86
  87 void
  88 metrics_postcalc(struct metrics *fs)
  89 {
  90   double sample_var;
  91   double cc = 0.0;
  92   double tc ;
  93   int k1, k2 ;
  94   int i;
  95   int j = 1;
  96
  97   struct weighted_value **data;
  98
  99
 100   int n_data;
 101
 102   fs->mean = fs->sum / fs->n;
 103
 104   sample_var = ( fs->ssq / fs->n  - fs->mean * fs->mean );
 105
 106   fs->var  = fs->n * sample_var / ( fs->n - 1) ;
 107   fs->stddev = sqrt(fs->var);
 108
 109
 110   /* FIXME: Check this is correct ???
 111      Shouldn't we use the sample variance ??? */
 112   fs->stderr = sqrt (fs->var / fs->n) ;
 113
 114   data = (struct weighted_value **) hsh_data(fs->ordered_data);
 115   n_data = hsh_count(fs->ordered_data);
 116
 117   fs->wv = xmalloc ( sizeof (struct weighted_value) * n_data);
 118
 119   for ( i = 0 ; i < n_data ; ++i )
 120     fs->wv[i] = *(data[i]);
 121
 122   sort (fs->wv, n_data, sizeof (struct weighted_value) ,
 123         (algo_compare_func *) compare_values, 0);
 124
 125
 126
 127   tc = fs->n * 0.05 ;
 128   k1 = -1;
 129   k2 = -1;
 130
 131
 132   for ( i = 0 ; i < n_data ; ++i )
 133     {
 134       cc += fs->wv[i].w;
 135       fs->wv[i].cc = cc;
 136
 137       fs->wv[i].rank = j + (fs->wv[i].w - 1) / 2.0 ;
 138
 139       j += fs->wv[i].w;
 140
 141       if ( cc < tc )
 142         k1 = i;
 143
 144     }
 145
 146   k2 = n_data;
 147   for ( i = n_data -1  ; i >= 0; --i )
 148     {
 149       if ( tc > fs->n - fs->wv[i].cc)
 150         k2 = i;
 151     }
 152
 153
 154   fs->trimmed_mean = 0;
 155   for ( i = k1 + 2 ; i <= k2 - 1 ; ++i )
 156     {
 157       fs->trimmed_mean += fs->wv[i].v.f * fs->wv[i].w;
 158     }
 159
 160
 161   fs->trimmed_mean += (fs->n - fs->wv[k2 - 1].cc - tc) * fs->wv[k2].v.f ;
 162   fs->trimmed_mean += (fs->wv[k1 + 1].cc - tc) * fs->wv[k1 + 1].v.f ;
 163   fs->trimmed_mean /= 0.9 * fs->n ;
 164
 165 }
 166
 167
 168 /* Functions for hashes */
 169
 170 void
 171 free_factor_stats(struct factor_statistics *f, int width UNUSED)
 172 {
 173   free (f);
 174 }
 175
 176 int
 177 compare_indep_values(const struct factor_statistics *f1,
 178                      const struct factor_statistics *f2,
 179                      int width)
 180 {
 181   return compare_values(f1->id, f2->id, width);
 182 }
 183
 184
 185 unsigned
 186 hash_indep_value(const struct factor_statistics *f, int width)
 187 {
 188   return hash_value(f->id, width);
 189 }