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