Set the filename on loading a new file
[pspp] / src / math / percentiles.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2008 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 "percentiles.h"
19 #include <math/order-stats.h>
20
21
22 #include "gettext.h"
23 #define _(msgid) gettext (msgid)
24 #define N_(msgid) msgid
25
26 #include <libpspp/assertion.h>
27 #include <data/val-type.h>
28 #include <gl/xalloc.h>
29 #include <data/variable.h>
30 #include <data/casereader.h>
31
32
33 const char *const ptile_alg_desc[] = {
34   "",
35   N_("HAverage"),
36   N_("Weighted Average"),
37   N_("Rounded"),
38   N_("Empirical"),
39   N_("Empirical with averaging")
40 };
41
42
43
44 double
45 percentile_calculate (const struct percentile *ptl, enum pc_alg alg)
46 {
47   struct percentile *mutable = (struct percentile *) ptl;
48   const struct order_stats *os = &ptl->parent;
49
50   assert (os->cc == ptl->w);
51
52   if ( ptl->g1 == SYSMIS)
53     mutable->g1 = (os->k[0].tc - os->k[0].cc) / os->k[0].c_p1;
54
55   if ( ptl->g1_star == SYSMIS)
56     mutable->g1_star = os->k[0].tc - os->k[0].cc;
57
58   if ( ptl->g2 == SYSMIS)
59     {
60       if ( os->k[1].c == 0 )
61         mutable->g2 = os->k[1].tc / os->k[1].c_p1;
62       else if ( os->k[1].c_p1 == 0 )
63         mutable->g2 = 0;
64       else
65         mutable->g2 = (os->k[1].tc - os->k[1].cc) / os->k[1].c_p1;
66     }
67
68   if ( ptl->g2_star == SYSMIS)
69     {
70       if ( os->k[1].c == 0 )
71         mutable->g2_star = os->k[1].tc;
72       else if ( os->k[1].c_p1 == 0 )
73         mutable->g2_star = 0;
74       else
75         mutable->g2_star = os->k[1].tc - os->k[1].cc;
76     }
77
78   switch (alg)
79     {
80     case PC_WAVERAGE:
81       if ( ptl->g1_star >= 1.0)
82         return os->k[0].y_p1;
83       else
84         {
85           double a = ( os->k[0].y == SYSMIS ) ? 0 : os->k[0].y;
86
87           if (os->k[0].c_p1 >= 1.0)
88             return (1 - ptl->g1_star) * a + ptl->g1_star * os->k[0].y_p1;
89           else
90             return (1 - ptl->g1) * a + ptl->g1 * os->k[0].y_p1;
91         }
92       break;
93
94     case PC_ROUND:
95       {
96         double a = ( os->k[0].y == SYSMIS ) ? 0 : os->k[0].y;
97
98         if (os->k[0].c_p1 >= 1.0)
99           return (ptl->g1_star < 0.5) ? a : os->k[0].y_p1;
100         else
101           return (ptl->g1 < 0.5) ? a : os->k[0].y_p1;
102       }
103       break;
104
105     case PC_EMPIRICAL:
106       if ( ptl->g1_star == 0 )
107         return os->k[0].y;
108       else
109         return os->k[0].y_p1;
110       break;
111
112     case PC_HAVERAGE:
113       if ( ptl->g2_star >= 1.0)
114         {
115           return os->k[1].y_p1;
116         }
117       else
118         {
119           double a = ( os->k[1].y == SYSMIS ) ? 0 : os->k[1].y;
120
121           if ( os->k[1].c_p1 >= 1.0)
122             {
123               if ( ptl->g2_star == 0)
124                 return os->k[1].y;
125
126               return (1 - ptl->g2_star) * a + ptl->g2_star * os->k[1].y_p1;
127             }
128           else
129             {
130               return (1 - ptl->g2) * a + ptl->g2 * os->k[1].y_p1;
131             }
132         }
133
134       break;
135
136     case PC_AEMPIRICAL:
137       if ( ptl->g1_star == 0 )
138         return (os->k[0].y + os->k[0].y_p1)/ 2.0;
139       else
140         return os->k[0].y_p1;
141       break;
142
143     default:
144       NOT_REACHED ();
145       break;
146     }
147
148   NOT_REACHED ();
149
150   return SYSMIS;
151 }
152
153
154 static void
155 destroy (struct statistic *stat)
156 {
157   struct order_stats *os = (struct order_stats *) stat;
158   free (os->k);
159   free (os);
160 }
161
162
163 struct order_stats *
164 percentile_create (double p, double W)
165 {
166   struct percentile *ptl = xzalloc (sizeof (*ptl));
167   struct order_stats *os = (struct order_stats *) ptl;
168   struct statistic *stat = (struct statistic *) ptl;
169
170   assert (p >= 0);
171   assert (p <= 1.0);
172
173   ptl->ptile = p;
174   ptl->w = W;
175
176   os->n_k = 2;
177   os->k = xcalloc (sizeof (*os->k), 2);
178   os->k[0].tc = W * p;
179   os->k[1].tc = (W + 1.0) * p;
180
181   ptl->g1 = ptl->g1_star = SYSMIS;
182   ptl->g2 = ptl->g2_star = SYSMIS;
183
184   os->k[1].y_p1 = os->k[1].y = SYSMIS;
185   os->k[0].y_p1 = os->k[0].y = SYSMIS;
186
187   stat->destroy = destroy;
188
189   return os;
190 }
191