1 /* PSPP - a program for statistical analysis.
2 Copyright (C) 2004 Free Software Foundation, Inc.
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.
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.
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/>. */
19 #include <libpspp/compiler.h>
20 #include "factor-stats.h"
21 #include "percentiles.h"
22 #include <libpspp/misc.h>
27 #define _(msgid) gettext (msgid)
28 #define N_(msgid) msgid
38 const char *const ptile_alg_desc[] = {
41 N_("Weighted Average"),
44 N_("Empirical with averaging")
50 /* Individual Percentile algorithms */
52 /* Closest observation to tc1 */
53 double ptile_round(const struct weighted_value **wv,
54 const struct ptile_params *par);
57 /* Weighted average at y_tc2 */
58 double ptile_haverage(const struct weighted_value **wv,
59 const struct ptile_params *par);
62 /* Weighted average at y_tc1 */
63 double ptile_waverage(const struct weighted_value **wv,
64 const struct ptile_params *par);
67 /* Empirical distribution function */
68 double ptile_empirical(const struct weighted_value **wv,
69 const struct ptile_params *par);
72 /* Empirical distribution function with averaging*/
73 double ptile_aempirical(const struct weighted_value **wv,
74 const struct ptile_params *par);
79 /* Closest observation to tc1 */
81 ptile_round(const struct weighted_value **wv,
82 const struct ptile_params *par)
90 if ( wv[par->k1 + 1]->w >= 1 )
92 if ( par->g1_star < 0.5 )
95 x = wv[par->k1 + 1]->v.f;
102 x = wv[par->k1 + 1]->v.f;
109 /* Weighted average at y_tc2 */
111 ptile_haverage(const struct weighted_value **wv,
112 const struct ptile_params *par)
117 if ( par->g2_star >= 1.0 )
118 return wv[par->k2 + 1]->v.f ;
120 /* Special case for k2 + 1 >= n_data
121 (actually it's not a special case, but just avoids indexing errors )
123 if ( par->g2_star == 0 )
125 assert(par->g2 == 0 );
126 return wv[par->k2]->v.f;
129 /* Ditto for k2 < 0 */
132 a = wv[par->k2]->v.f;
135 if ( wv[par->k2 + 1]->w >= 1.0 )
136 return ( (1 - par->g2_star) * a +
137 par->g2_star * wv[par->k2 + 1]->v.f);
139 return ( (1 - par->g2) * a +
140 par->g2 * wv[par->k2 + 1]->v.f);
146 /* Weighted average at y_tc1 */
148 ptile_waverage(const struct weighted_value **wv,
149 const struct ptile_params *par)
153 if ( par->g1_star >= 1.0 )
154 return wv[par->k1 + 1]->v.f ;
158 a = wv[par->k1]->v.f;
161 if ( wv[par->k1 + 1]->w >= 1.0 )
162 return ( (1 - par->g1_star) * a +
163 par->g1_star * wv[par->k1 + 1]->v.f);
165 return ( (1 - par->g1) * a +
166 par->g1 * wv[par->k1 + 1]->v.f);
170 /* Empirical distribution function */
172 ptile_empirical(const struct weighted_value **wv,
173 const struct ptile_params *par)
175 if ( par->g1_star > 0 )
176 return wv[par->k1 + 1]->v.f;
178 return wv[par->k1]->v.f;
183 /* Empirical distribution function with averageing */
185 ptile_aempirical(const struct weighted_value **wv,
186 const struct ptile_params *par)
188 if ( par->g1_star > 0 )
189 return wv[par->k1 + 1]->v.f;
191 return (wv[par->k1]->v.f + wv[par->k1 + 1]->v.f ) / 2.0 ;
196 /* Compute the percentile p */
197 double ptile(double p,
198 const struct weighted_value **wv,
201 enum pc_alg algorithm);
207 const struct weighted_value **wv,
210 enum pc_alg algorithm)
216 struct ptile_params pp;
226 for ( i = 0 ; i < n_data ; ++i )
228 if ( wv[i]->cc <= tc1 )
231 if ( wv[i]->cc <= tc2 )
239 pp.g1 = ( tc1 - wv[pp.k1]->cc ) / wv[pp.k1 + 1]->w;
240 pp.g1_star = tc1 - wv[pp.k1]->cc ;
244 pp.g1 = tc1 / wv[pp.k1 + 1]->w;
249 if ( pp.k2 + 1 >= n_data )
258 pp.g2 = ( tc2 - wv[pp.k2]->cc ) / wv[pp.k2 + 1]->w;
259 pp.g2_star = tc2 - wv[pp.k2]->cc ;
263 pp.g2 = tc2 / wv[pp.k2 + 1]->w;
271 result = ptile_haverage(wv, &pp);
274 result = ptile_waverage(wv, &pp);
277 result = ptile_round(wv, &pp);
280 result = ptile_empirical(wv, &pp);
283 result = ptile_aempirical(wv, &pp);
294 Calculate the values of the percentiles in pc_hash.
295 wv is a sorted array of weighted values of the data set.
298 ptiles(struct hsh_table *pc_hash,
299 const struct weighted_value **wv,
302 enum pc_alg algorithm)
304 struct hsh_iterator hi;
305 struct percentile *p;
309 for ( p = hsh_first(pc_hash, &hi);
311 p = hsh_next(pc_hash, &hi))
313 p->v = ptile(p->p/100.0 , wv, n_data, w, algorithm);
319 /* Calculate Tukey's Hinges */
321 tukey_hinges(const struct weighted_value **wv,
328 double c_star = DBL_MAX;
334 for ( i = 0 ; i < n_data ; ++i )
336 c_star = MIN(c_star, wv[i]->w);
339 if ( c_star > 1 ) c_star = 1;
341 d = floor((w/c_star + 3 ) / 2.0)/ 2.0;
344 l[1] = w/2.0 + c_star/2.0;
345 l[2] = w + c_star - d*c_star;
351 for ( i = 0 ; i < n_data ; ++i )
353 if ( l[0] >= wv[i]->cc ) h[0] = i ;
354 if ( l[1] >= wv[i]->cc ) h[1] = i ;
355 if ( l[2] >= wv[i]->cc ) h[2] = i ;
358 for ( i = 0 ; i < 3 ; i++ )
362 a_star = l[i] - wv[h[i]]->cc ;
366 if ( h[i] + 1 >= n_data )
368 assert( a_star < 1 ) ;
369 hinge[i] = (1 - a_star) * wv[h[i]]->v.f;
374 a = a_star / ( wv[h[i] + 1]->cc ) ;
379 hinge[i] = wv[h[i] + 1]->v.f ;
383 if ( wv[h[i] + 1]->w >= 1)
385 hinge[i] = ( 1 - a_star) * wv[h[i]]->v.f
386 + a_star * wv[h[i] + 1]->v.f;
391 hinge[i] = (1 - a) * wv[h[i]]->v.f + a * wv[h[i] + 1]->v.f;
395 assert(hinge[0] <= hinge[1]);
396 assert(hinge[1] <= hinge[2]);
402 ptile_compare(const struct percentile *p1,
403 const struct percentile *p2,
411 else if (p1->p < p2->p)
420 ptile_hash(const struct percentile *p, void *aux UNUSED)
422 return hsh_hash_double(p->p);