src/math/percentiles.c

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