gui: Allow File|Open to select an encoding for system files.
[pspp] / 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_data_rw_idx (cp, NP_IDX_Y)->f = y;
67   case_data_rw_idx (cp, NP_IDX_NS)->f = ns;
68   case_data_rw_idx (cp, NP_IDX_DNS)->f = dns;
69   casewriter_write (np->writer, cp);
70
71   np->prev_cc = cc;
72 }
73
74 struct np *
75 np_create (double n, double mean, double var)
76 {
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   np->n = n;
86   np->mean = mean;
87
88   np->stddev = sqrt (var);
89
90   np->y_min = np->ns_min = np->dns_min = DBL_MAX;
91   np->y_max = np->ns_max = np->dns_max = -DBL_MAX;
92
93   proto = caseproto_create ();
94   for (i = 0; i < n_NP_IDX; i++)
95     proto = caseproto_add_width (proto, 0);
96   np->writer = autopaging_writer_create (proto);
97   caseproto_unref (proto);
98
99   os->k = 0;
100   stat->destroy = destroy;
101   stat->accumulate = acc;
102
103   return np;
104 }