1 /* PSPP - A program for statistical analysis . -*-c-*-
3 Copyright (C) 2004 Free Software Foundation, Inc.
4 Author: John Darrington 2004
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.
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.
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
22 #include "factor_stats.h"
23 #include "percentiles.h"
27 #define _(msgid) gettext (msgid)
28 #define N_(msgid) msgid
41 const char *ptile_alg_desc[] = {
44 N_("Weighted Average"),
47 N_("Empirical with averaging")
53 /* Individual Percentile algorithms */
55 /* Closest observation to tc1 */
56 double ptile_round(const struct weighted_value **wv,
57 const struct ptile_params *par);
60 /* Weighted average at y_tc2 */
61 double ptile_haverage(const struct weighted_value **wv,
62 const struct ptile_params *par);
65 /* Weighted average at y_tc1 */
66 double ptile_waverage(const struct weighted_value **wv,
67 const struct ptile_params *par);
70 /* Empirical distribution function */
71 double ptile_empirical(const struct weighted_value **wv,
72 const struct ptile_params *par);
75 /* Empirical distribution function with averaging*/
76 double ptile_aempirical(const struct weighted_value **wv,
77 const struct ptile_params *par);
82 /* Closest observation to tc1 */
84 ptile_round(const struct weighted_value **wv,
85 const struct ptile_params *par)
93 if ( wv[par->k1 + 1]->w >= 1 )
95 if ( par->g1_star < 0.5 )
98 x = wv[par->k1 + 1]->v.f;
105 x = wv[par->k1 + 1]->v.f;
112 /* Weighted average at y_tc2 */
114 ptile_haverage(const struct weighted_value **wv,
115 const struct ptile_params *par)
120 if ( par->g2_star >= 1.0 )
121 return wv[par->k2 + 1]->v.f ;
123 /* Special case for k2 + 1 >= n_data
124 (actually it's not a special case, but just avoids indexing errors )
126 if ( par->g2_star == 0 )
128 assert(par->g2 == 0 );
129 return wv[par->k2]->v.f;
132 /* Ditto for k2 < 0 */
135 a = wv[par->k2]->v.f;
138 if ( wv[par->k2 + 1]->w >= 1.0 )
139 return ( (1 - par->g2_star) * a +
140 par->g2_star * wv[par->k2 + 1]->v.f);
142 return ( (1 - par->g2) * a +
143 par->g2 * wv[par->k2 + 1]->v.f);
149 /* Weighted average at y_tc1 */
151 ptile_waverage(const struct weighted_value **wv,
152 const struct ptile_params *par)
156 if ( par->g1_star >= 1.0 )
157 return wv[par->k1 + 1]->v.f ;
161 a = wv[par->k1]->v.f;
164 if ( wv[par->k1 + 1]->w >= 1.0 )
165 return ( (1 - par->g1_star) * a +
166 par->g1_star * wv[par->k1 + 1]->v.f);
168 return ( (1 - par->g1) * a +
169 par->g1 * wv[par->k1 + 1]->v.f);
173 /* Empirical distribution function */
175 ptile_empirical(const struct weighted_value **wv,
176 const struct ptile_params *par)
178 if ( par->g1_star > 0 )
179 return wv[par->k1 + 1]->v.f;
181 return wv[par->k1]->v.f;
186 /* Empirical distribution function with averageing */
188 ptile_aempirical(const struct weighted_value **wv,
189 const struct ptile_params *par)
191 if ( par->g1_star > 0 )
192 return wv[par->k1 + 1]->v.f;
194 return (wv[par->k1]->v.f + wv[par->k1 + 1]->v.f ) / 2.0 ;
199 /* Compute the percentile p */
200 double ptile(double p,
201 const struct weighted_value **wv,
204 enum pc_alg algorithm);
210 const struct weighted_value **wv,
213 enum pc_alg algorithm)
219 struct ptile_params pp;
229 for ( i = 0 ; i < n_data ; ++i )
231 if ( wv[i]->cc <= tc1 )
234 if ( wv[i]->cc <= tc2 )
242 pp.g1 = ( tc1 - wv[pp.k1]->cc ) / wv[pp.k1 + 1]->w;
243 pp.g1_star = tc1 - wv[pp.k1]->cc ;
247 pp.g1 = tc1 / wv[pp.k1 + 1]->w;
252 if ( pp.k2 + 1 >= n_data )
261 pp.g2 = ( tc2 - wv[pp.k2]->cc ) / wv[pp.k2 + 1]->w;
262 pp.g2_star = tc2 - wv[pp.k2]->cc ;
266 pp.g2 = tc2 / wv[pp.k2 + 1]->w;
274 result = ptile_haverage(wv, &pp);
277 result = ptile_waverage(wv, &pp);
280 result = ptile_round(wv, &pp);
283 result = ptile_empirical(wv, &pp);
286 result = ptile_aempirical(wv, &pp);
297 Calculate the values of the percentiles in pc_hash.
298 wv is a sorted array of weighted values of the data set.
301 ptiles(struct hsh_table *pc_hash,
302 const struct weighted_value **wv,
305 enum pc_alg algorithm)
307 struct hsh_iterator hi;
308 struct percentile *p;
312 for ( p = hsh_first(pc_hash, &hi);
314 p = hsh_next(pc_hash, &hi))
316 p->v = ptile(p->p/100.0 , wv, n_data, w, algorithm);
322 /* Calculate Tukey's Hinges */
324 tukey_hinges(const struct weighted_value **wv,
331 double c_star = DBL_MAX;
337 for ( i = 0 ; i < n_data ; ++i )
339 c_star = min(c_star, wv[i]->w);
342 if ( c_star > 1 ) c_star = 1;
344 d = floor((w/c_star + 3 ) / 2.0)/ 2.0;
347 l[1] = w/2.0 + c_star/2.0;
348 l[2] = w + c_star - d*c_star;
354 for ( i = 0 ; i < n_data ; ++i )
356 if ( l[0] >= wv[i]->cc ) h[0] = i ;
357 if ( l[1] >= wv[i]->cc ) h[1] = i ;
358 if ( l[2] >= wv[i]->cc ) h[2] = i ;
361 for ( i = 0 ; i < 3 ; i++ )
365 a_star = l[i] - wv[h[i]]->cc ;
369 if ( h[i] + 1 >= n_data )
371 assert( a_star < 1 ) ;
372 hinge[i] = (1 - a_star) * wv[h[i]]->v.f;
377 a = a_star / ( wv[h[i] + 1]->cc ) ;
382 hinge[i] = wv[h[i] + 1]->v.f ;
386 if ( wv[h[i] + 1]->w >= 1)
388 hinge[i] = ( 1 - a_star) * wv[h[i]]->v.f
389 + a_star * wv[h[i] + 1]->v.f;
394 hinge[i] = (1 - a) * wv[h[i]]->v.f + a * wv[h[i] + 1]->v.f;
398 assert(hinge[0] <= hinge[1]);
399 assert(hinge[1] <= hinge[2]);
405 ptile_compare(const struct percentile *p1,
406 const struct percentile *p2,
414 else if (p1->p < p2->p)
423 ptile_hash(const struct percentile *p, void *aux UNUSED)
425 return hsh_hash_double(p->p);