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