pintos-os.org Git - pspp/blob - src/math/trimmed-mean.c

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