Fixed bug #22073
[pspp] / src / math / percentiles.c
1 /* PSPP - a program for statistical analysis.
2    Copyright (C) 2004 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 <assert.h>
19 #include <data/val-type.h>
20 #include <libpspp/compiler.h>
21 #include "factor-stats.h"
22 #include "percentiles.h"
23 #include <libpspp/misc.h>
24
25 #include "minmax.h"
26
27 #include "gettext.h"
28 #define _(msgid) gettext (msgid)
29 #define N_(msgid) msgid
30
31 struct ptile_params
32 {
33   double g1, g1_star;
34   double g2, g2_star;
35   int k1, k2;
36 };
37
38
39 const char *const ptile_alg_desc[] = {
40   "",
41   N_("HAverage"),
42   N_("Weighted Average"),
43   N_("Rounded"),
44   N_("Empirical"),
45   N_("Empirical with averaging")
46 };
47
48
49
50
51 /* Individual Percentile algorithms */
52
53 /* Closest observation to tc1 */
54 double ptile_round(const struct weighted_value **wv,
55                    const struct ptile_params *par);
56
57
58 /* Weighted average at y_tc2 */
59 double ptile_haverage(const struct weighted_value **wv,
60                    const struct ptile_params *par);
61
62
63 /* Weighted average at y_tc1 */
64 double ptile_waverage(const struct weighted_value **wv,
65                    const struct ptile_params *par);
66
67
68 /* Empirical distribution function */
69 double ptile_empirical(const struct weighted_value **wv,
70                    const struct ptile_params *par);
71
72
73 /* Empirical distribution function with averaging*/
74 double ptile_aempirical(const struct weighted_value **wv,
75                    const struct ptile_params *par);
76
77
78
79
80 /* Closest observation to tc1 */
81 double
82 ptile_round(const struct weighted_value **wv,
83             const struct ptile_params *par)
84 {
85   double x;
86   double a=0;
87
88   if ( par->k1 >= 0 )
89     a = wv[par->k1]->v.f;
90
91   if ( wv[par->k1 + 1]->w >= 1 )
92     {
93       if ( par->g1_star < 0.5 )
94         x = a;
95       else
96         x = wv[par->k1 + 1]->v.f;
97     }
98   else
99     {
100       if ( par->g1 < 0.5 )
101         x = a;
102       else
103         x = wv[par->k1 + 1]->v.f;
104
105     }
106
107   return x;
108 }
109
110 /* Weighted average at y_tc2 */
111 double
112 ptile_haverage(const struct weighted_value **wv,
113                const struct ptile_params *par)
114 {
115
116   double a=0;
117
118   if ( par->g2_star >= 1.0 )
119       return wv[par->k2 + 1]->v.f ;
120
121   /* Special case  for k2 + 1 >= n_data
122      (actually it's not a special case, but just avoids indexing errors )
123    */
124   if ( par->g2_star == 0 )
125     {
126       assert(par->g2 == 0 );
127       return wv[par->k2]->v.f;
128     }
129
130   /* Ditto for k2 < 0 */
131   if ( par->k2 >= 0 )
132     {
133       a = wv[par->k2]->v.f;
134     }
135
136   if ( wv[par->k2 + 1]->w >= 1.0 )
137     return ( (1 - par->g2_star) *  a   +
138              par->g2_star * wv[par->k2 + 1]->v.f);
139   else
140     return ( (1 - par->g2) * a +
141              par->g2 * wv[par->k2 + 1]->v.f);
142
143 }
144
145
146
147 /* Weighted average at y_tc1 */
148 double
149 ptile_waverage(const struct weighted_value **wv,
150                const struct ptile_params *par)
151 {
152   double a=0;
153
154   if ( par->g1_star >= 1.0 )
155       return wv[par->k1 + 1]->v.f ;
156
157   if ( par->k1 >= 0 )
158     {
159       a = wv[par->k1]->v.f;
160     }
161
162   if ( wv[par->k1 + 1]->w >= 1.0 )
163     return ( (1 - par->g1_star) * a +
164              par->g1_star * wv[par->k1 + 1]->v.f);
165   else
166     return ( (1 - par->g1) * a +
167              par->g1 * wv[par->k1 + 1]->v.f);
168 }
169
170
171 /* Empirical distribution function */
172 double
173 ptile_empirical(const struct weighted_value **wv,
174                const struct ptile_params *par)
175 {
176   if ( par->g1_star > 0 )
177     return wv[par->k1 + 1]->v.f;
178   else
179     return wv[par->k1]->v.f;
180 }
181
182
183
184 /* Empirical distribution function with averageing */
185 double
186 ptile_aempirical(const struct weighted_value **wv,
187                const struct ptile_params *par)
188 {
189   if ( par->g1_star > 0 )
190     return wv[par->k1 + 1]->v.f;
191   else
192     return (wv[par->k1]->v.f + wv[par->k1 + 1]->v.f ) / 2.0 ;
193 }
194
195
196
197 /* Compute the percentile p */
198 double ptile(double p,
199              const struct weighted_value **wv,
200              int n_data,
201              double w,
202              enum pc_alg algorithm);
203
204
205
206 double
207 ptile(double p,
208       const struct weighted_value **wv,
209       int n_data,
210       double w,
211       enum pc_alg algorithm)
212 {
213   int i;
214   double tc1, tc2;
215   double result;
216
217   struct ptile_params pp;
218
219   assert( p <= 1.0);
220
221   tc1 = w * p ;
222   tc2 = (w + 1) * p ;
223
224   pp.k1 = -1;
225   pp.k2 = -1;
226
227   for ( i = 0 ; i < n_data ; ++i )
228     {
229       if ( wv[i]->cc <= tc1 )
230         pp.k1 = i;
231
232       if ( wv[i]->cc <= tc2 )
233         pp.k2 = i;
234
235     }
236
237
238   if ( pp.k1 >= 0 )
239     {
240       pp.g1 = ( tc1 - wv[pp.k1]->cc ) / wv[pp.k1 + 1]->w;
241       pp.g1_star = tc1 -  wv[pp.k1]->cc ;
242     }
243   else
244     {
245       pp.g1 = tc1 / wv[pp.k1 + 1]->w;
246       pp.g1_star = tc1 ;
247     }
248
249
250   if ( pp.k2  + 1 >= n_data )
251     {
252       pp.g2 = 0 ;
253       pp.g2_star = 0;
254     }
255   else
256     {
257       if ( pp.k2 >= 0 )
258         {
259           pp.g2 = ( tc2 - wv[pp.k2]->cc ) / wv[pp.k2 + 1]->w;
260           pp.g2_star = tc2 -  wv[pp.k2]->cc ;
261         }
262       else
263         {
264           pp.g2 = tc2 / wv[pp.k2 + 1]->w;
265           pp.g2_star = tc2 ;
266         }
267     }
268
269   switch ( algorithm )
270     {
271     case PC_HAVERAGE:
272       result = ptile_haverage(wv, &pp);
273       break;
274     case PC_WAVERAGE:
275       result = ptile_waverage(wv, &pp);
276       break;
277     case PC_ROUND:
278       result = ptile_round(wv, &pp);
279       break;
280     case PC_EMPIRICAL:
281       result = ptile_empirical(wv, &pp);
282       break;
283     case PC_AEMPIRICAL:
284       result = ptile_aempirical(wv, &pp);
285       break;
286     default:
287       result = SYSMIS;
288     }
289
290   return result;
291 }
292
293
294 /*
295    Calculate the values of the percentiles in pc_hash.
296    wv is  a sorted array of weighted values of the data set.
297 */
298 void
299 ptiles(struct hsh_table *pc_hash,
300        const struct weighted_value **wv,
301        int n_data,
302        double w,
303        enum pc_alg algorithm)
304 {
305   struct hsh_iterator hi;
306   struct percentile *p;
307
308   if ( !pc_hash )
309     return ;
310   for ( p = hsh_first(pc_hash, &hi);
311         p != 0 ;
312         p = hsh_next(pc_hash, &hi))
313     {
314       p->v = ptile(p->p/100.0 , wv, n_data, w, algorithm);
315     }
316
317 }
318
319
320 /* Calculate Tukey's Hinges */
321 void
322 tukey_hinges(const struct weighted_value **wv,
323              int n_data,
324              double w,
325              double hinge[3]
326              )
327 {
328   int i;
329   double c_star = DBL_MAX;
330   double d;
331   double l[3];
332   int h[3];
333   double a, a_star;
334
335   for ( i = 0 ; i < n_data ; ++i )
336     {
337       c_star = MIN(c_star, wv[i]->w);
338     }
339
340   if ( c_star > 1 ) c_star = 1;
341
342   d = floor((w/c_star + 3 ) / 2.0)/ 2.0;
343
344   l[0] = d*c_star;
345   l[1] = w/2.0 + c_star/2.0;
346   l[2] = w + c_star - d*c_star;
347
348   h[0]=-1;
349   h[1]=-1;
350   h[2]=-1;
351
352   for ( i = 0 ; i < n_data ; ++i )
353     {
354       if ( l[0] >= wv[i]->cc ) h[0] = i ;
355       if ( l[1] >= wv[i]->cc ) h[1] = i ;
356       if ( l[2] >= wv[i]->cc ) h[2] = i ;
357     }
358
359   for ( i = 0 ; i < 3 ; i++ )
360     {
361
362       if ( h[i] >= 0 )
363         a_star = l[i] - wv[h[i]]->cc ;
364       else
365         a_star = l[i];
366
367       if ( h[i] + 1 >= n_data )
368       {
369               assert( a_star < 1 ) ;
370               hinge[i] = (1 - a_star) * wv[h[i]]->v.f;
371               continue;
372       }
373       else
374       {
375               a = a_star / ( wv[h[i] + 1]->cc ) ;
376       }
377
378       if ( a_star >= 1.0 )
379         {
380           hinge[i] = wv[h[i] + 1]->v.f ;
381           continue;
382         }
383
384       if ( wv[h[i] + 1]->w >= 1)
385         {
386           hinge[i] = ( 1 - a_star) * wv[h[i]]->v.f
387             + a_star * wv[h[i] + 1]->v.f;
388
389           continue;
390         }
391
392       hinge[i] = (1 - a) * wv[h[i]]->v.f + a * wv[h[i] + 1]->v.f;
393
394     }
395
396   assert(hinge[0] <= hinge[1]);
397   assert(hinge[1] <= hinge[2]);
398
399 }
400
401
402 int
403 ptile_compare(const struct percentile *p1,
404                    const struct percentile *p2,
405                    void *aux UNUSED)
406 {
407
408   int cmp;
409
410   if ( p1->p == p2->p)
411     cmp = 0 ;
412   else if (p1->p < p2->p)
413     cmp = -1 ;
414   else
415     cmp = +1;
416
417   return cmp;
418 }
419
420 unsigned
421 ptile_hash(const struct percentile *p, void *aux UNUSED)
422 {
423   return hsh_hash_double(p->p);
424 }
425
426