Oops. Neglected to add new files.

[pspp-builds.git] / src / math / trimmed-mean.c

diff --git a/src/math/trimmed-mean.c b/src/math/trimmed-mean.c
new file mode 100644 (file)
index 0000000..da3d424
--- /dev/null
+++ b/src/math/trimmed-mean.c
@@ -0,0 +1,94 @@
+/* PSPP - a program for statistical analysis.
+   Copyright (C) 2008 Free Software Foundation, Inc.
+
+   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 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, see <http://www.gnu.org/licenses/>. */
+
+#include <config.h>
+#include "trimmed-mean.h"
+#include <math/order-stats.h>
+
+#include <gl/xalloc.h>
+#include <libpspp/assertion.h>
+#include <math.h>
+#include <data/val-type.h>
+
+
+static void
+acc (struct statistic *s, const struct ccase *cx UNUSED, double c, double cc, double y)
+{
+  struct trimmed_mean *tm = (struct trimmed_mean *) s;
+  struct order_stats *os = (struct order_stats *) s;
+
+  if ( cc > os->k[0].tc && cc < os->k[1].tc)
+      tm->sum += c * y;
+
+  if ( tm->cyk1p1 == SYSMIS && cc >os->k[0].tc)
+      tm->cyk1p1 = c * y;
+}
+
+static void
+destroy (struct statistic *s)
+{
+  struct order_stats *os = (struct order_stats *) s;
+  free (os->k);
+  free (s);
+}
+
+struct statistic *
+trimmed_mean_create (double W, double tail)
+{
+  struct trimmed_mean *tm = xzalloc (sizeof (*tm));
+  struct order_stats *os = (struct order_stats *) tm;
+  struct statistic *stat = (struct statistic *) tm;
+
+  os->n_k = 2;
+  os->k = xcalloc (sizeof (*os->k), 2);
+
+  assert (tail >= 0);
+  assert (tail <= 1);
+
+  os->k[0].tc = tail * W;
+  os->k[1].tc = W * (1 - tail);
+
+  stat->accumulate = acc;
+  stat->destroy = destroy;
+
+  tm->cyk1p1 = SYSMIS;
+  tm->w = W;
+  tm->tail = tail;
+
+  return stat;
+}
+
+
+double
+trimmed_mean_calculate (const struct trimmed_mean *tm)
+{
+  const struct order_stats *os = (const struct order_stats *) tm;
+
+  assert (os->cc == tm->w);
+
+  return
+    (
+     (os->k[0].cc_p1 - os->k[0].tc) * os->k[0].y_p1
+     -
+     (os->k[1].cc - os->k[1].tc) * os->k[1].y_p1
+     +
+     tm->sum
+     -
+     tm->cyk1p1
+     )
+    /
+    ( (1.0 - 2 * tm->tail) * tm->w);
+}