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