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