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