Oops. Neglected to add new files.
authorJohn Darrington <jmd@pc-188.(none)>
Fri, 5 Sep 2008 02:04:40 +0000 (10:04 +0800)
committerJohn Darrington <jmd@pc-188.(none)>
Fri, 5 Sep 2008 02:04:40 +0000 (10:04 +0800)
Actually I'm *sure* that I did do "git add" on all these,
but for whatever reason git didn't check them in.
Hopefully it's right this time.  One day I'll work out
how to use git effectively.  One day ...

14 files changed:
src/math/box-whisker.c [new file with mode: 0644]
src/math/box-whisker.h [new file with mode: 0644]
src/math/extrema.c [new file with mode: 0644]
src/math/extrema.h [new file with mode: 0644]
src/math/np.c [new file with mode: 0644]
src/math/np.h [new file with mode: 0644]
src/math/order-stats.c [new file with mode: 0644]
src/math/order-stats.h [new file with mode: 0644]
src/math/statistic.h [new file with mode: 0644]
src/math/trimmed-mean.c [new file with mode: 0644]
src/math/trimmed-mean.h [new file with mode: 0644]
src/math/tukey-hinges.c [new file with mode: 0644]
src/math/tukey-hinges.h [new file with mode: 0644]
tests/bugs/examine-crash.sh [new file with mode: 0755]

diff --git a/src/math/box-whisker.c b/src/math/box-whisker.c
new file mode 100644 (file)
index 0000000..288fc07
--- /dev/null
@@ -0,0 +1,139 @@
+/* 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 "box-whisker.h"
+#include "order-stats.h"
+#include "tukey-hinges.h"
+#include <gl/xalloc.h>
+#include <libpspp/assertion.h>
+#include <math.h>
+#include <float.h>
+#include <data/val-type.h>
+#include <libpspp/str.h>
+#include <data/case.h>
+#include <data/variable.h>
+
+static void
+destroy (struct statistic *s)
+{
+  struct order_stats *os = (struct order_stats *) s;
+  struct box_whisker *bw = (struct box_whisker *) s;
+  struct ll *ll;
+
+  for (ll = ll_head (&bw->outliers); ll != ll_null (&bw->outliers); )
+    {
+      struct outlier *e = ll_data (ll, struct outlier, ll);
+
+      ll = ll_next (ll);
+
+      ds_destroy (&e->label);
+      free (e);
+    }
+
+  free (os->k);
+  free (s);
+};
+
+
+static void
+acc (struct statistic *s, const struct ccase *cx,
+     double c UNUSED, double cc UNUSED, double y)
+{
+  struct box_whisker *bw = (struct box_whisker *) s;
+  bool extreme;
+  struct outlier *o;
+
+  if ( y < bw->hinges[2] + bw->step)
+      bw->whiskers[1] = y;
+
+  if (bw->whiskers[0] == SYSMIS ||  bw->hinges[0] - bw->step > y)
+      bw->whiskers[0] = y;
+
+  if ( y > bw->hinges[2] + bw->step)
+    extreme = (y > bw->hinges[2] + 2 * bw->step) ;
+
+  else if (y < bw->hinges[0] - bw->step)
+    extreme = (y < bw->hinges[0] - 2 * bw->step) ;
+
+  else
+    return;
+
+  o = xzalloc (sizeof *o) ;
+  o->value = y;
+  o->extreme = extreme;
+  ds_init_empty (&o->label);
+
+  if (bw->id_var)
+    var_append_value_name (bw->id_var,
+                          case_data (cx, bw->id_var),
+                          &o->label);
+  else
+    ds_put_format (&o->label,
+                  "%ld",
+                  (casenumber) case_data_idx (cx, bw->casenumber_idx)->f);
+
+  ll_push_head (&bw->outliers, &o->ll);
+}
+
+void
+box_whisker_whiskers (const struct box_whisker *bw, double whiskers[2])
+{
+  whiskers[0] = bw->whiskers[0];
+  whiskers[1] = bw->whiskers[1];
+}
+
+void
+box_whisker_hinges (const struct box_whisker *bw, double hinges[3])
+{
+  hinges[0] = bw->hinges[0];
+  hinges[1] = bw->hinges[1];
+  hinges[2] = bw->hinges[2];
+}
+
+const struct ll_list *
+box_whisker_outliers (const struct box_whisker *bw)
+{
+  return &bw->outliers;
+}
+
+struct statistic *
+box_whisker_create (const struct tukey_hinges *th,
+                   const struct variable *id_var,  size_t casenumber_idx)
+{
+  struct box_whisker *w = xzalloc (sizeof (*w));
+  struct order_stats *os = (struct order_stats *) w;
+  struct statistic *stat = (struct statistic *) w;
+
+  os->n_k = 0;
+
+  stat->destroy = destroy;
+  stat->accumulate = acc;
+
+  tukey_hinges_calculate (th, w->hinges);
+
+  w->casenumber_idx = casenumber_idx;
+  w->id_var = id_var;
+
+  w->step = (w->hinges[2] - w->hinges[0]) * 1.5;
+
+  w->whiskers[1] = w->hinges[2];
+  w->whiskers[0] = SYSMIS;
+
+  ll_init (&w->outliers);
+
+  return stat;
+}
diff --git a/src/math/box-whisker.h b/src/math/box-whisker.h
new file mode 100644 (file)
index 0000000..5202b64
--- /dev/null
@@ -0,0 +1,65 @@
+/* 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/>. */
+
+#ifndef __MATH_BOX_WHISKER_H__
+#define __MATH_BOX_WHISKER_H__
+
+#include <stddef.h>
+#include <libpspp/ll.h>
+#include <libpspp/str.h>
+#include "order-stats.h"
+
+/* This module calculates the statistics typically displayed by box-plots.
+   However, there's no reason not to use it for other purposes too.
+ */
+struct tukey_hinges;
+
+
+struct outlier
+{
+  double value;
+  struct string label;
+  bool extreme;
+  struct ll ll;
+};
+
+
+struct box_whisker
+{
+  struct order_stats parent;
+
+  double hinges[3];
+  double whiskers[2];
+
+  struct ll_list outliers;
+
+  double step;
+
+  size_t casenumber_idx;
+  const struct variable *id_var;
+};
+
+struct statistic * box_whisker_create (const struct tukey_hinges *,
+                                        const struct variable *, size_t);
+
+void box_whisker_whiskers (const struct box_whisker *bw, double whiskers[2]);
+
+void box_whisker_hinges (const struct box_whisker *bw, double hinges[2]);
+
+const struct ll_list * box_whisker_outliers (const struct box_whisker *bw);
+
+
+#endif
diff --git a/src/math/extrema.c b/src/math/extrema.c
new file mode 100644 (file)
index 0000000..617c7ac
--- /dev/null
@@ -0,0 +1,144 @@
+/* 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 "extrema.h"
+#include <xalloc.h>
+#include <data/case.h>
+#include <data/val-type.h>
+#include <libpspp/compiler.h>
+#include <libpspp/ll.h>
+#include <stdlib.h>
+
+struct extrema
+{
+  size_t capacity;
+  size_t n;
+  struct ll_list list;
+
+  ll_compare_func *cmp_func;
+};
+
+
+static int
+cmp_descending (const struct ll *a_, const struct ll *b_, void *aux UNUSED)
+{
+  const struct extremum *a = ll_data (a_, struct extremum, ll);
+  const struct extremum *b = ll_data (b_, struct extremum, ll);
+
+  if ( a->value > b->value) return -1;
+
+  return (a->value < b->value);
+}
+
+static int
+cmp_ascending (const struct ll *a_, const struct ll *b_, void *aux UNUSED)
+{
+  const struct extremum *a = ll_data (a_, struct extremum, ll);
+  const struct extremum *b = ll_data (b_, struct extremum, ll);
+
+  if ( a->value < b->value) return -1;
+
+  return (a->value > b->value);
+}
+
+
+struct extrema *
+extrema_create (size_t n, enum extreme_end end)
+{
+  struct extrema *extrema = xzalloc (sizeof *extrema);
+  extrema->capacity = n;
+
+  if ( end == EXTREME_MAXIMA )
+    extrema->cmp_func = cmp_descending;
+  else
+    extrema->cmp_func = cmp_ascending;
+
+  ll_init (&extrema->list);
+
+  return extrema;
+}
+
+void
+extrema_destroy (struct extrema *extrema)
+{
+  struct ll *ll = ll_head (&extrema->list);
+
+  while (ll != ll_null (&extrema->list))
+    {
+      struct extremum *e = ll_data (ll, struct extremum, ll);
+
+      ll = ll_next (ll);
+      free (e);
+    }
+
+  free (extrema);
+}
+
+
+void
+extrema_add (struct extrema *extrema, double val,
+            double weight,
+            casenumber location)
+{
+  struct extremum *e = xzalloc (sizeof *e) ;
+  e->value = val;
+  e->location = location;
+  e->weight = weight;
+
+  if ( val == SYSMIS)
+    {
+      free (e);
+      return;
+    }
+
+  ll_insert_ordered (ll_head (&extrema->list), ll_null (&extrema->list),
+                      &e->ll,  extrema->cmp_func, NULL);
+
+  if ( extrema->n++ > extrema->capacity)
+    {
+      struct ll *tail = ll_tail (&extrema->list);
+      struct extremum *et = ll_data (tail, struct extremum, ll);
+
+      ll_remove (tail);
+
+      free (et);
+    }
+}
+
+
+const struct ll_list *
+extrema_list (const struct extrema *ex)
+{
+  return &ex->list;
+}
+
+
+bool
+extrema_top (const struct extrema *ex, double *v)
+{
+  const struct extremum  *top;
+
+  if ( ll_is_empty (&ex->list))
+    return false;
+
+  top = (const struct extremum *)
+    ll_data (ll_head(&ex->list), struct extremum, ll);
+
+  *v = top->value;
+
+  return true;
+}
diff --git a/src/math/extrema.h b/src/math/extrema.h
new file mode 100644 (file)
index 0000000..d891c53
--- /dev/null
@@ -0,0 +1,58 @@
+/* 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/>. */
+
+#ifndef __EXTREMA_H__
+#define __EXTREMA_H__ 1
+
+#include <stddef.h>
+#include <data/case.h>
+#include <libpspp/ll.h>
+
+struct extremum
+{
+  double value;
+  casenumber location;
+  double weight;
+
+  /* Internal use only */
+  struct ll ll;
+};
+
+
+enum extreme_end
+  {
+    EXTREME_MAXIMA,
+    EXTREME_MINIMA
+  };
+
+struct extrema;
+
+struct extrema *extrema_create (size_t n, enum extreme_end);
+
+void extrema_destroy (struct extrema *extrema);
+
+void extrema_add (struct extrema *extrema, double val,
+                 double weight,
+                 casenumber location);
+
+void extrema_show (const struct extrema *extrema);
+
+const struct ll_list * extrema_list (const struct extrema *);
+
+bool extrema_top (const struct extrema *, double *);
+
+
+#endif
diff --git a/src/math/np.c b/src/math/np.c
new file mode 100644 (file)
index 0000000..e189b47
--- /dev/null
@@ -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 "np.h"
+#include <math/moments.h>
+#include <gl/xalloc.h>
+#include <stdlib.h>
+#include <math.h>
+#include <gsl/gsl_cdf.h>
+#include <libpspp/compiler.h>
+#include <data/case.h>
+#include <data/casewriter.h>
+
+static void
+destroy (struct statistic *stat)
+{
+  struct order_stats *os = (struct order_stats *) stat;
+  free (os);
+}
+
+
+static void
+acc (struct statistic *s, const struct ccase *cx UNUSED,
+     double c, double cc, double y)
+{
+  struct ccase cp;
+  struct np *np = (struct np *) s;
+  double rank = np->prev_cc + (c + 1) / 2.0;
+
+  double ns = gsl_cdf_ugaussian_Pinv (rank / ( np->n + 1 ));
+
+  double z = (y - np->mean) / np->stddev;
+
+  double dns = z - ns;
+
+  maximize (&np->ns_max, ns);
+  minimize (&np->ns_min, ns);
+
+  maximize (&np->dns_max, dns);
+  minimize (&np->dns_min, dns);
+
+  maximize (&np->y_max, y);
+  minimize (&np->y_min, y);
+
+  case_create (&cp, n_NP_IDX);
+
+  case_data_rw_idx (&cp, NP_IDX_Y)->f = y;
+  case_data_rw_idx (&cp, NP_IDX_NS)->f = ns;
+  case_data_rw_idx (&cp, NP_IDX_DNS)->f = dns;
+
+  casewriter_write (np->writer, &cp);
+
+  np->prev_cc = cc;
+}
+
+struct order_stats *
+np_create (const struct moments1 *m)
+{
+  double variance;
+  struct np *np = xzalloc (sizeof (*np));
+  struct statistic *stat = (struct statistic *) np;
+  struct order_stats *os = (struct order_stats *) np;
+
+  np->prev_cc = 0;
+
+  moments1_calculate (m, &np->n, &np->mean, &variance, NULL, NULL);
+
+  np->stddev = sqrt (variance);
+
+  np->y_min = np->ns_min = np->dns_min = DBL_MAX;
+  np->y_max = np->ns_max = np->dns_max = -DBL_MAX;
+
+  np->writer = autopaging_writer_create (n_NP_IDX);
+
+  os->k = 0;
+  stat->destroy = destroy;
+  stat->accumulate = acc;
+
+  return os;
+}
diff --git a/src/math/np.h b/src/math/np.h
new file mode 100644 (file)
index 0000000..7db51f7
--- /dev/null
@@ -0,0 +1,59 @@
+/* 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/>. */
+
+#ifndef __NP_H__
+#define __NP_H__
+
+#include "order-stats.h"
+
+struct moments1;
+struct casewriter;
+
+enum
+  {
+    NP_IDX_Y = 0,
+    NP_IDX_NS,
+    NP_IDX_DNS,
+    n_NP_IDX
+  };
+
+struct np
+{
+  struct order_stats parent;
+
+  double n;
+  double mean;
+  double stddev;
+
+
+  double prev_cc;
+
+  double ns_min;
+  double ns_max;
+
+  double dns_min;
+  double dns_max;
+
+  double y_min;
+  double y_max;
+
+  struct casewriter *writer;
+};
+
+
+struct order_stats * np_create (const struct moments1 *);
+
+#endif
diff --git a/src/math/order-stats.c b/src/math/order-stats.c
new file mode 100644 (file)
index 0000000..ca4160f
--- /dev/null
@@ -0,0 +1,159 @@
+/* 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 "order-stats.h"
+#include <libpspp/assertion.h>
+#include <data/val-type.h>
+#include <gl/xalloc.h>
+#include <data/variable.h>
+#include <data/casereader.h>
+#include <string.h>
+
+#if 0
+
+#include <stdio.h>
+
+static void
+order_stats_dump_k1 (const struct order_stats *os)
+{
+  struct k *k = &os->k[0];
+  printf ("K1: tc %g; c %g cc %g ccp %g\n",
+         k->tc, k->c, k->cc, k->cc_p1);
+
+}
+
+static void
+order_stats_dump_k2 (const struct order_stats *os)
+{
+  struct k *k = &os->k[1];
+  printf ("K2: tc %g; c %g cc %g ccp %g\n",
+         k->tc, k->c, k->cc, k->cc_p1);
+}
+
+
+void
+order_stats_dump (const struct order_stats *os)
+{
+  order_stats_dump_k1 (os);
+  order_stats_dump_k2 (os);
+}
+
+#endif
+
+static void
+update_k_lower (struct k *kk,
+               double y_i, double c_i, double cc_i)
+{
+  if ( cc_i <= kk->tc )
+    {
+      kk->cc = cc_i;
+      kk->c = c_i;
+      kk->y = y_i;
+    }
+}
+
+
+static void
+update_k_upper (struct k *kk,
+               double y_i, double c_i, double cc_i)
+{
+  if ( cc_i > kk->tc && kk->c_p1 == 0)
+    {
+      kk->cc_p1 = cc_i;
+      kk->c_p1 = c_i;
+      kk->y_p1 = y_i;
+    }
+}
+
+
+static void
+update_k_values (const struct ccase *cx, double y_i, double c_i, double cc_i,
+                struct order_stats **os, size_t n_os)
+{
+  int j;
+
+  for (j = 0 ; j < n_os ; ++j)
+    {
+      int k;
+      struct order_stats *tos = os[j];
+      struct statistic  *stat = (struct statistic *) tos;
+      for (k = 0 ; k < tos->n_k; ++k)
+       {
+         struct k *myk = &tos->k[k];
+         update_k_lower (myk, y_i, c_i, cc_i);
+         update_k_upper (myk, y_i, c_i, cc_i);
+       }
+
+      if ( stat->accumulate )
+       stat->accumulate (stat, cx, c_i, cc_i, y_i);
+
+      tos->cc = cc_i;
+    }
+}
+
+
+void
+order_stats_accumulate (struct order_stats **os, size_t nos,
+                       struct casereader *reader,
+                       const struct variable *wv,
+                       const struct variable *var,
+                       enum mv_class exclude)
+{
+  struct ccase cx;
+  struct ccase prev_cx;
+  double prev_value = -DBL_MAX;
+
+  double cc_i = 0;
+  double c_i = 0;
+
+  case_nullify (&prev_cx);
+
+  for (; casereader_read (reader, &cx); case_destroy (&cx))
+    {
+      const double weight = wv ? case_data (&cx, wv)->f : 1.0;
+      const double this_value = case_data (&cx, var)->f;
+
+      /* The casereader MUST be sorted */
+      assert (this_value >= prev_value);
+
+      if ( var_is_value_missing (var, case_data (&cx, var), exclude))
+       continue;
+
+      case_destroy (&prev_cx);
+
+      if ( prev_value == -DBL_MAX || prev_value == this_value)
+       c_i += weight;
+
+      if ( prev_value > -DBL_MAX && this_value > prev_value)
+       {
+         update_k_values (&prev_cx, prev_value, c_i, cc_i, os, nos);
+         c_i = weight;
+       }
+
+      cc_i += weight;
+      prev_value = this_value;
+      case_clone (&prev_cx, &cx);
+    }
+
+  update_k_values (&prev_cx, prev_value, c_i, cc_i, os, nos);
+  case_destroy (&prev_cx);
+
+  casereader_destroy (reader);
+}
+
+
+
diff --git a/src/math/order-stats.h b/src/math/order-stats.h
new file mode 100644 (file)
index 0000000..cea50ed
--- /dev/null
@@ -0,0 +1,60 @@
+/* PSPP - a program for statistical analysis.
+   Copyright (C) 2004, 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/>. */
+
+#ifndef __ORDER_STATS_H__
+#define __ORDER_STATS_H__
+
+#include <stddef.h>
+#include <math/statistic.h>
+
+struct casereader;
+struct variable;
+
+/*
+  cc <= tc < cc_p1
+*/
+struct k
+{
+  double tc;
+  double cc;
+  double cc_p1;
+  double c;
+  double c_p1;
+  double y;
+  double y_p1;
+};
+
+
+struct order_stats
+{
+  struct statistic parent;
+  int n_k;
+  struct k *k;
+
+  double cc;
+};
+
+enum mv_class;
+
+void order_stats_dump (const struct order_stats *os);
+
+void order_stats_accumulate (struct order_stats **ptl, size_t nos,
+                            struct casereader *reader,
+                            const struct variable *wv,
+                            const struct variable *var,
+                            enum mv_class exclude);
+
+#endif
diff --git a/src/math/statistic.h b/src/math/statistic.h
new file mode 100644 (file)
index 0000000..987264b
--- /dev/null
@@ -0,0 +1,40 @@
+/* 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/>. */
+
+#ifndef __STATISTIC_H__
+#define __STATISTIC_H__
+
+#include <stddef.h>
+
+struct ccase ;
+
+struct statistic
+{
+  void (*accumulate) (struct statistic *, const struct ccase *cx, double c, double cc, double y);
+  void (*destroy) (struct statistic *);
+};
+
+static inline void statistic_destroy (struct statistic *s);
+
+
+static inline void
+statistic_destroy (struct statistic *s)
+{
+  if (s) s->destroy (s);
+}
+
+
+#endif
diff --git a/src/math/trimmed-mean.c b/src/math/trimmed-mean.c
new file mode 100644 (file)
index 0000000..da3d424
--- /dev/null
@@ -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);
+}
diff --git a/src/math/trimmed-mean.h b/src/math/trimmed-mean.h
new file mode 100644 (file)
index 0000000..9339cab
--- /dev/null
@@ -0,0 +1,42 @@
+/* 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/>. */
+
+#ifndef __TRIMMED_MEAN_H__
+#define __TRIMMED_MEAN_H__
+
+#include <stddef.h>
+
+#include "order-stats.h"
+
+
+struct trimmed_mean
+{
+  struct order_stats parent;
+
+  /* The partial sum */
+  double sum;
+
+  /* The product of c_{k_1+1} and y_{k_1 + 1} */
+  double cyk1p1;
+
+  double w;
+  double tail;
+};
+
+struct statistic * trimmed_mean_create (double W, double c_min);
+double trimmed_mean_calculate (const struct trimmed_mean *);
+
+#endif
diff --git a/src/math/tukey-hinges.c b/src/math/tukey-hinges.c
new file mode 100644 (file)
index 0000000..95a79c1
--- /dev/null
@@ -0,0 +1,101 @@
+/* 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 "tukey-hinges.h"
+#include <math/order-stats.h>
+
+#include <gl/xalloc.h>
+#include <libpspp/assertion.h>
+#include <math.h>
+
+void
+tukey_hinges_calculate (const struct tukey_hinges *th, double hinge[3])
+{
+  double a[3];
+  double a_star[3];
+  int i;
+  const struct order_stats *os = &th->parent;
+
+  for (i = 0 ; i < 3 ; ++i)
+    {
+      a_star[i] = os->k[i].tc - os->k[i].cc;
+      a[i] = a_star[i] / os->k[i].c_p1;
+
+      if (a_star[i] < 1)
+       {
+         if (os->k[i].c_p1 >= 1 )
+           {
+             hinge[i] = (1 - a_star[i]) * os->k[i].y
+               + a_star[i] * os->k[i].y_p1;
+           }
+         else
+           {
+             hinge[i] = (1 - a[i]) * os->k[i].y
+               + a[i] * os->k[i].y_p1;
+           }
+       }
+      else
+       {
+         hinge[i] = os->k[i].y_p1;
+       }
+
+    }
+}
+
+static void
+destroy (struct statistic *s)
+{
+  struct order_stats *os = (struct order_stats *) s;
+
+  free (os->k);
+  free (s);
+};
+
+struct statistic *
+tukey_hinges_create (double W, double c_min)
+{
+  double d;
+  struct tukey_hinges *th = xzalloc (sizeof (*th));
+  struct order_stats *os = (struct order_stats *) th;
+  struct statistic *stat = (struct statistic *) th;
+
+  assert (c_min >= 0);
+
+  os->n_k = 3;
+  os->k = xcalloc (sizeof (*os->k), 3);
+
+  if ( c_min >= 1.0)
+    {
+      d = floor ((W + 3) / 2.0) / 2.0;
+
+      os->k[0].tc = d;
+      os->k[1].tc = W/2.0 + 0.5;
+      os->k[2].tc = W + 1 - d;
+    }
+  else
+    {
+      d = floor ((W/c_min + 3.0)/ 2.0) / 2.0 ;
+      os->k[0].tc = d * c_min;
+      os->k[1].tc = (W + c_min) / 2.0;
+      os->k[2].tc = W + c_min * (1 - d);
+    }
+
+
+  stat->destroy = destroy;
+
+  return stat;
+}
diff --git a/src/math/tukey-hinges.h b/src/math/tukey-hinges.h
new file mode 100644 (file)
index 0000000..d87691f
--- /dev/null
@@ -0,0 +1,37 @@
+/* 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/>. */
+
+#ifndef __TUKEY_HINGES_H__
+#define __TUKEY_HINGES_H__
+
+#include <stddef.h>
+
+#include "order-stats.h"
+
+
+struct tukey_hinges
+{
+  struct order_stats parent;
+};
+
+struct statistic * tukey_hinges_create (double W, double c_min);
+
+
+void tukey_hinges_calculate (const struct tukey_hinges *h, double hinge[3]);
+
+
+
+#endif
diff --git a/tests/bugs/examine-crash.sh b/tests/bugs/examine-crash.sh
new file mode 100755 (executable)
index 0000000..6cd172f
--- /dev/null
@@ -0,0 +1,80 @@
+#!/bin/sh
+
+# This program tests for a bug which crashed EXAMINE
+
+TEMPDIR=/tmp/pspp-tst-$$
+TESTFILE=$TEMPDIR/`basename $0`.sps
+
+# ensure that top_srcdir and top_builddir  are absolute
+if [ -z "$top_srcdir" ] ; then top_srcdir=. ; fi
+if [ -z "$top_builddir" ] ; then top_builddir=. ; fi
+top_srcdir=`cd $top_srcdir; pwd`
+top_builddir=`cd $top_builddir; pwd`
+
+PSPP=$top_builddir/src/ui/terminal/pspp
+
+STAT_CONFIG_PATH=$top_srcdir/config
+export STAT_CONFIG_PATH
+
+LANG=C
+export LANG
+
+
+cleanup()
+{
+     if [ x"$PSPP_TEST_NO_CLEANUP" != x ] ; then 
+       echo "NOT cleaning $TEMPDIR"
+       return ; 
+     fi
+     rm -rf $TEMPDIR
+}
+
+
+fail()
+{
+    echo $activity
+    echo FAILED
+    cleanup;
+    exit 1;
+}
+
+
+no_result()
+{
+    echo $activity
+    echo NO RESULT;
+    cleanup;
+    exit 2;
+}
+
+pass()
+{
+    cleanup;
+    exit 0;
+}
+
+mkdir -p $TEMPDIR
+
+cd $TEMPDIR
+
+cat <<EOF > $TESTFILE
+data list list /a * x * y *.
+begin data.
+3 1 3
+5 1 4
+7 2 3
+end data.
+
+examine a by x by y
+       /statistics=DESCRIPTIVES
+       . 
+EOF
+if [ $? -ne 0 ] ; then no_result ; fi
+
+
+activity="run program"
+$SUPERVISOR $PSPP --testing-mode -o raw-ascii $TESTFILE
+if [ $? -ne 0 ] ; then fail ; fi
+
+
+pass;