src/math/percentiles.c

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