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