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