pintos-os.org Git - pspp/blob - src/math/np.c

   1 /* PSPP - a program for statistical analysis.
   2    Copyright (C) 2008, 2009, 2011 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
  19 #include "math/np.h"
  20
  21 #include <gsl/gsl_cdf.h>
  22 #include <math.h>
  23 #include <stdlib.h>
  24
  25 #include "data/case.h"
  26 #include "data/casewriter.h"
  27 #include "libpspp/cast.h"
  28 #include "libpspp/compiler.h"
  29 #include "libpspp/misc.h"
  30 #include "math/moments.h"
  31
  32 #include "gl/xalloc.h"
  33
  34 static void
  35 destroy (struct statistic *stat)
  36 {
  37   struct np *np = UP_CAST (stat, struct np, parent.parent);
  38   free (np);
  39 }
  40
  41
  42 static void
  43 acc (struct statistic *s, const struct ccase *cx UNUSED,
  44      double c, double cc, double y)
  45 {
  46   struct ccase *cp;
  47   struct np *np = UP_CAST (s, struct np, parent.parent);
  48   double rank = np->prev_cc + (c + 1) / 2.0;
  49
  50   double ns = gsl_cdf_ugaussian_Pinv (rank / (np->n + 1));
  51
  52   double z = (y - np->mean) / np->stddev;
  53
  54   double dns = z - ns;
  55
  56   maximize (&np->ns_max, ns);
  57   minimize (&np->ns_min, ns);
  58
  59   maximize (&np->dns_max, dns);
  60   minimize (&np->dns_min, dns);
  61
  62   maximize (&np->y_max, y);
  63   minimize (&np->y_min, y);
  64
  65   cp = case_create (casewriter_get_proto (np->writer));
  66   *case_num_rw_idx (cp, NP_IDX_Y) = y;
  67   *case_num_rw_idx (cp, NP_IDX_NS) = ns;
  68   *case_num_rw_idx (cp, NP_IDX_DNS) = dns;
  69   casewriter_write (np->writer, cp);
  70
  71   np->prev_cc = cc;
  72 }
  73
  74 /* Creates and returns a data structure whose accumulated results can be used
  75    to produce a normal probability plot.  The caller must supply the weighted
  76    sample size N and the mean MEAN and variance VAR of the distribution, then
  77    feed in the data with order_stats_accumulate() or
  78    order_stats_accumulate_idx().
  79
  80    There is no function to produce the results, which appear in "struct np" for
  81    passing directly to np_plot_create() or dnp_plot_create().
  82
  83    The caller must eventually destroy the returned structure, with
  84    statistic_destroy(). */
  85 struct np *
  86 np_create (double n, double mean, double var)
  87 {
  88   struct np *np = XZALLOC (struct np);
  89   struct order_stats *os = &np->parent;
  90   struct statistic *stat = &os->parent;
  91   struct caseproto *proto;
  92   int i;
  93
  94   np->prev_cc = 0;
  95
  96   np->n = n;
  97   np->mean = mean;
  98
  99   np->stddev = sqrt (var);
 100
 101   np->y_min = np->ns_min = np->dns_min = DBL_MAX;
 102   np->y_max = np->ns_max = np->dns_max = -DBL_MAX;
 103
 104   proto = caseproto_create ();
 105   for (i = 0; i < n_NP_IDX; i++)
 106     proto = caseproto_add_width (proto, 0);
 107   np->writer = autopaging_writer_create (proto);
 108   caseproto_unref (proto);
 109
 110   os->k = 0;
 111   stat->destroy = destroy;
 112   stat->accumulate = acc;
 113
 114   return np;
 115 }