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