src/math/trimmed-mean.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2008, 2009 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 #include "trimmed-mean.h"
  19 #include <math/order-stats.h>
  20
  21 #include <gl/xalloc.h>
  22 #include <libpspp/assertion.h>
  23 #include <libpspp/cast.h>
  24 #include <math.h>
  25 #include <data/val-type.h>
  26
  27
  28 static void
  29 acc (struct statistic *s, const struct ccase *cx UNUSED, double c, double cc, double y)
  30 {
  31   struct trimmed_mean *tm = UP_CAST (s, struct trimmed_mean, parent.parent);
  32   struct order_stats *os = &tm->parent;
  33
  34   if ( cc > os->k[0].tc && cc < os->k[1].tc)
  35       tm->sum += c * y;
  36
  37   if ( tm->cyk1p1 == SYSMIS && cc >os->k[0].tc)
  38       tm->cyk1p1 = c * y;
  39 }
  40
  41 static void
  42 destroy (struct statistic *s)
  43 {
  44   struct trimmed_mean *tm = UP_CAST (s, struct trimmed_mean, parent.parent);
  45   struct order_stats *os = &tm->parent;
  46   free (os->k);
  47   free (tm);
  48 }
  49
  50 struct trimmed_mean *
  51 trimmed_mean_create (double W, double tail)
  52 {
  53   struct trimmed_mean *tm = xzalloc (sizeof (*tm));
  54   struct order_stats *os = &tm->parent;
  55   struct statistic *stat = &os->parent;
  56
  57   os->n_k = 2;
  58   os->k = xcalloc (sizeof (*os->k), 2);
  59
  60   assert (tail >= 0);
  61   assert (tail <= 1);
  62
  63   os->k[0].tc = tail * W;
  64   os->k[1].tc = W * (1 - tail);
  65
  66   stat->accumulate = acc;
  67   stat->destroy = destroy;
  68
  69   tm->cyk1p1 = SYSMIS;
  70   tm->w = W;
  71   tm->tail = tail;
  72
  73   return tm;
  74 }
  75
  76
  77 double
  78 trimmed_mean_calculate (const struct trimmed_mean *tm)
  79 {
  80   const struct order_stats *os = (const struct order_stats *) tm;
  81
  82   assert (os->cc == tm->w);
  83
  84   return
  85     (
  86      (os->k[0].cc_p1 - os->k[0].tc) * os->k[0].y_p1
  87      -
  88      (os->k[1].cc - os->k[1].tc) * os->k[1].y_p1
  89      +
  90      tm->sum
  91      -
  92      tm->cyk1p1
  93      )
  94     /
  95     ( (1.0 - 2 * tm->tail) * tm->w);
  96 }