src/math/np.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 "np.h"
  19
  20 #include <gsl/gsl_cdf.h>
  21 #include <math.h>
  22 #include <stdlib.h>
  23
  24 #include <data/case.h>
  25 #include <data/casewriter.h>
  26 #include <libpspp/compiler.h>
  27 #include <libpspp/cast.h>
  28 #include <libpspp/misc.h>
  29 #include <math/moments.h>
  30
  31 #include "xalloc.h"
  32
  33 static void
  34 destroy (struct statistic *stat)
  35 {
  36   struct np *np = UP_CAST (stat, struct np, parent.parent);
  37   free (np);
  38 }
  39
  40
  41 static void
  42 acc (struct statistic *s, const struct ccase *cx UNUSED,
  43      double c, double cc, double y)
  44 {
  45   struct ccase *cp;
  46   struct np *np = UP_CAST (s, struct np, parent.parent);
  47   double rank = np->prev_cc + (c + 1) / 2.0;
  48
  49   double ns = gsl_cdf_ugaussian_Pinv (rank / ( np->n + 1 ));
  50
  51   double z = (y - np->mean) / np->stddev;
  52
  53   double dns = z - ns;
  54
  55   maximize (&np->ns_max, ns);
  56   minimize (&np->ns_min, ns);
  57
  58   maximize (&np->dns_max, dns);
  59   minimize (&np->dns_min, dns);
  60
  61   maximize (&np->y_max, y);
  62   minimize (&np->y_min, y);
  63
  64   cp = case_create (casewriter_get_proto (np->writer));
  65   case_data_rw_idx (cp, NP_IDX_Y)->f = y;
  66   case_data_rw_idx (cp, NP_IDX_NS)->f = ns;
  67   case_data_rw_idx (cp, NP_IDX_DNS)->f = dns;
  68   casewriter_write (np->writer, cp);
  69
  70   np->prev_cc = cc;
  71 }
  72
  73 struct np *
  74 np_create (const struct moments1 *m)
  75 {
  76   double variance;
  77   struct np *np = xzalloc (sizeof (*np));
  78   struct order_stats *os = &np->parent;
  79   struct statistic *stat = &os->parent;
  80   struct caseproto *proto;
  81   int i;
  82
  83   np->prev_cc = 0;
  84
  85   moments1_calculate (m, &np->n, &np->mean, &variance, NULL, NULL);
  86
  87   np->stddev = sqrt (variance);
  88
  89   np->y_min = np->ns_min = np->dns_min = DBL_MAX;
  90   np->y_max = np->ns_max = np->dns_max = -DBL_MAX;
  91
  92   proto = caseproto_create ();
  93   for (i = 0; i < n_NP_IDX; i++)
  94     proto = caseproto_add_width (proto, 0);
  95   np->writer = autopaging_writer_create (proto);
  96   caseproto_unref (proto);
  97
  98   os->k = 0;
  99   stat->destroy = destroy;
 100   stat->accumulate = acc;
 101
 102   return np;
 103 }