--- /dev/null
+/* 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;
+}
--- /dev/null
+/* 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
--- /dev/null
+/* 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;
+}
--- /dev/null
+/* 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
--- /dev/null
+/* 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;
+}
--- /dev/null
+/* 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
--- /dev/null
+/* 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);
+}
+
+
+
--- /dev/null
+/* 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
--- /dev/null
+/* 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
--- /dev/null
+/* 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);
+}
--- /dev/null
+/* 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
--- /dev/null
+/* 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;
+}
--- /dev/null
+/* 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
--- /dev/null
+#!/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;