c_cpp/etc/calc/cal/statistics.cal

   1 /*
   2  *  statistics - Some assorted statistics functions.
   3  *
   4  * Copyright (C) 2013 Christoph Zurnieden
   5  *
   6  * statistics is open software; you can redistribute it and/or modify it under
   7  * the terms of the version 2.1 of the GNU Lesser General Public License
   8  * as published by the Free Software Foundation.
   9  *
  10  * statistics is distributed in the hope that it will be useful, but WITHOUT
  11  * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
  12  * or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU Lesser General
  13  * Public License for more details.
  14  *
  15  * A copy of version 2.1 of the GNU Lesser General Public License is
  16  * distributed with calc under the filename COPYING-LGPL.  You should have
  17  * received a copy with calc; if not, write to Free Software Foundation, Inc.
  18  * 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.
  19  *
  20  * @(#) $Revision: 30.3 $
  21  * @(#) $Id: statistics.cal,v 30.3 2013/08/11 08:41:38 chongo Exp $
  22  * @(#) $Source: /usr/local/src/bin/calc/cal/RCS/statistics.cal,v $
  23  *
  24  * Under source code control:   2013/08/11 01:31:28
  25  * File existed as early as:    2013
  26  */
  27
  28
  29 static resource_debug_level;
  30 resource_debug_level = config("resource_debug", 0);
  31
  32
  33 /*
  34   get dependencies
  35 */
  36 read -once factorial2 brentsolve
  37
  38
  39 /*******************************************************************************
  40  *
  41  *
  42  *                Continuous distributions
  43  *
  44  *
  45  ******************************************************************************/
  46
  47 /* regularized incomplete gamma function like in Octave, hence the name */
  48 define gammaincoctave(z,a){
  49   local tmp;
  50   tmp = gamma(z);
  51   return (tmp-gammainc(a,z))/tmp;
  52 }
  53
  54 /* Inverse incomplete beta function. Old and slow. */
  55 static __CZ__invbeta_a;
  56 static __CZ__invbeta_b;
  57 static __CZ__invbeta_x;
  58 define __CZ__invbeta(x){
  59    return __CZ__invbeta_x-__CZ__ibetaas63(x,__CZ__invbeta_a,__CZ__invbeta_b);
  60 }
  61
  62 define invbetainc_slow(x,a,b){
  63   local flag ret eps;
  64   /* place checks and balances here */
  65   eps = epsilon();
  66   if(.5 < x){
  67     __CZ__invbeta_x = 1 - x;
  68     __CZ__invbeta_a = b;
  69     __CZ__invbeta_b = a;
  70     flag = 1;
  71   }
  72   else{
  73     __CZ__invbeta_x = x;
  74     __CZ__invbeta_a = a;
  75     __CZ__invbeta_b = b;
  76     flag = 0;
  77   }
  78
  79   ret =  brentsolve2(0,1,1);
  80
  81   if(flag == 1)
  82     ret = 1-ret;
  83   epsilon(eps);
  84   return ret;
  85 }
  86
  87 /* Inverse incomplete beta function. Still old but not as slow as the function
  88    above. */
  89 /*
  90   Purpose:
  91
  92     invbetainc computes inverse of the incomplete Beta function.
  93
  94   Licensing:
  95
  96     This code is distributed under the GNU LGPL license.
  97
  98   Modified:
  99
 100     10 August 2013
 101
 102   Author:
 103
 104     Original FORTRAN77 version by GW Cran, KJ Martin, GE Thomas.
 105     C version by John Burkardt.
 106     Calc version by Christoph Zurnieden
 107
 108   Reference:
 109
 110     GW Cran, KJ Martin, GE Thomas,
 111     Remark AS R19 and Algorithm AS 109:
 112     A Remark on Algorithms AS 63: The Incomplete Beta Integral
 113     and AS 64: Inverse of the Incomplete Beta Integeral,
 114     Applied Statistics,
 115     Volume 26, Number 1, 1977, pages 111-114.
 116
 117   Parameters:
 118
 119     Input,  P, Q, the parameters of the incomplete
 120     Beta function.
 121
 122     Input,  BETA, the logarithm of the value of
 123     the complete Beta function.
 124
 125     Input,  ALPHA, the value of the incomplete Beta
 126     function.  0 <= ALPHA <= 1.
 127
 128     Output,  the argument of the incomplete
 129     Beta function which produces the value ALPHA.
 130
 131   Local Parameters:
 132
 133     Local, SAE, the most negative decimal exponent
 134     which does not cause an underflow.
 135 */
 136 define invbetainc(x,a,b){
 137   return __CZ__invbetainc(a,b,lnbeta(a,b),x);
 138 }
 139
 140 define __CZ__invbetainc(p,q,beta,alpha){
 141   local a acu adj fpu g h iex indx pp prev qq r s sae sq t tx value;
 142   local w xin y yprev places eps;
 143
 144   /* Dirty trick, don't try at home */
 145   eps= epsilon(epsilon()^2);
 146   sae = -((log(1/epsilon())/log(2))//2);
 147   fpu = 10.0^sae;
 148
 149   places = highbit(1 + int(1/epsilon())) + 1;
 150   value = alpha;
 151   if( p <= 0.0 ){
 152     epsilon(eps);
 153     return newerror("invbeta: argument p <= 0");
 154   }
 155   if( q <= 0.0 ){
 156     epsilon(eps);
 157     return newerror("invbeta: argument q <= 0");
 158   }
 159
 160   if( alpha < 0.0 || 1.0 < alpha ){
 161     epsilon(eps);
 162     return newerror("invbeta: argument alpha out of domain");
 163   }
 164   if( alpha == 0.0 ){
 165     epsilon(eps);
 166     return 0;
 167   }
 168   if( alpha == 1.0 ){
 169     epsilon(eps);
 170     return 1;
 171   }
 172   if ( 0.5 < alpha ){
 173     a = 1.0 - alpha;
 174     pp = q;
 175     qq = p;
 176     indx = 1;
 177   }
 178   else{
 179     a = alpha;
 180     pp = p;
 181     qq = q;
 182     indx = 0;
 183   }
 184   r = sqrt ( - ln ( a * a ) );
 185
 186   y = r-(2.30753+0.27061*r)/(1.0+(0.99229+0.04481*r)*r);
 187
 188   if ( 1.0 < pp && 1.0 < qq ){
 189     r = ( y * y - 3.0 ) / 6.0;
 190     s = 1.0 / ( pp + pp - 1.0 );
 191     t = 1.0 / ( qq + qq - 1.0 );
 192     h = 2.0 / ( s + t );
 193     w = y*sqrt(h+r)/h-(t-s)*(r+5.0/6.0-2.0/(3.0*h));
 194     value = pp / ( pp + qq * exp ( w + w ) );
 195   }
 196   else{
 197     r = qq + qq;
 198     t = 1.0 / ( 9.0 * qq );
 199     t = r *  ( 1.0 - t + y * sqrt ( t )^ 3 );
 200
 201     if ( t <= 0.0 ){
 202       value = 1.0 - exp ( ( ln ( ( 1.0 - a ) * qq ) + beta ) / qq );
 203     }
 204     else{
 205       t = ( 4.0 * pp + r - 2.0 ) / t;
 206
 207       if ( t <= 1.0 ) {
 208         value = exp ( ( ln ( a * pp ) + beta ) / pp );
 209       }
 210       else{
 211         value = 1.0 - 2.0 / ( t + 1.0 );
 212       }
 213     }
 214   }
 215   r = 1.0 - pp;
 216   t = 1.0 - qq;
 217   yprev = 0.0;
 218   sq = 1.0;
 219   prev = 1.0;
 220
 221   if ( value < 0.0001 )
 222     value = 0.0001;
 223
 224   if ( 0.9999 < value )
 225     value = 0.9999;
 226
 227   acu = 10^sae;
 228
 229   for ( ; ; ){
 230     y = bround(__CZ__ibetaas63( value, pp, qq, beta),places);
 231     xin = value;
 232     y = bround(exp(ln(y-a)+(beta+r*ln(xin)+t*ln(1.0- xin ) )),places);
 233
 234     if ( y * yprev <= 0.0 ) {
 235       prev = max ( sq, fpu );
 236     }
 237
 238     g = 1.0;
 239
 240     for ( ; ; ){
 241       for ( ; ; ){
 242         adj = g * y;
 243         sq = adj * adj;
 244         if ( sq < prev ){
 245           tx = value - adj;
 246           if ( 0.0 <= tx && tx <= 1.0 )  break;
 247         }
 248         g = g / 3.0;
 249       }
 250       if ( prev <= acu ){
 251         if ( indx )
 252           value = 1.0 - value;
 253         epsilon(eps);
 254         return value;
 255       }
 256       if ( y * y <= acu ){
 257         if ( indx )
 258           value = 1.0 - value;
 259         epsilon(eps);
 260         return value;
 261       }
 262       if ( tx != 0.0 && tx != 1.0 )
 263         break;
 264       g = g / 3.0;
 265     }
 266     if ( tx == value )  break;
 267     value = tx;
 268     yprev = y;
 269   }
 270   if ( indx )
 271     value = 1.0 - value;
 272
 273   epsilon(eps);
 274   return value;
 275 }
 276
 277 /*******************************************************************************
 278  *
 279  *
 280  *                Beta distribution
 281  *
 282  *
 283  ******************************************************************************/
 284
 285 define betapdf(x,a,b){
 286   if(x<0 || x>1) return newerror("betapdf: parameter x out of domain");
 287   if(a<=0) return newerror("betapdf: parameter a out of domain");
 288   if(b<=0) return newerror("betapdf: parameter b out of domain");
 289
 290   return 1/beta(a,b) *x^(a-1)*(1-x)^(b-1);
 291 }
 292
 293 define betacdf(x,a,b){
 294   if(x<0 || x>1) return newerror("betacdf: parameter x out of domain");
 295   if(a<=0) return newerror("betacdf: parameter a out of domain");
 296   if(b<=0) return newerror("betacdf: parameter b out of domain");
 297
 298   return betainc(x,a,b);
 299 }
 300
 301 define betacdfinv(x,a,b){
 302   return invbetainc(x,a,b);
 303 }
 304
 305 define betamedian(a,b){
 306   local t106 t104 t103 t105 approx ret;
 307   if(a == b) return 1/2;
 308   if(a == 1 && b > 0) return 1-(1/2)^(1/b);
 309   if(a > 0 && b == 1) return   (1/2)^(1/a);
 310   if(a == 3 && b == 2){
 311     /* Yes, the author is not ashamed to ask Maxima for the exact solution
 312        of a quartic equation. */
 313     t103 = ( (2^(3/2))/27 +4/27  )^(1/3);
 314     t104 = sqrt( ( 9*t103^2 + 4*t103 + 2  )/(t103)  )/3;
 315     t105 = -t103-2/(9*t103) +8/9;
 316     t106 = sqrt( (27*t104*t105+16)/(t104)  )/(2*3^(3/2));
 317     return -t106+t104/2+1/3;
 318   }
 319   if(a == 2 && b == 3){
 320     t103 = ( (2^(3/2))/27 +4/27  )^(1/3);
 321     t104 = sqrt( ( 9*t103^2 + 4*t103 + 2  )/(t103)  )/3;
 322     t105 = -t103-2/(9*t103) +8/9;
 323     t106 = sqrt( (27*t104*t105+16)/(t104)  )/(2*3^(3/2));
 324     return 1-(-t106+t104/2+1/3);
 325   }
 326   return invbetainc(1/2,a,b);
 327 }
 328
 329 define betamode(a,b){
 330   if(a + b == 2) return newerror("betamod: a + b = 2 = division by zero");
 331   return (a-1)/(a+b-2);
 332 }
 333
 334 define betavariance(a,b){
 335   return (a*b)/( (a+b)^2*(a+b+1) );
 336 }
 337
 338 define betalnvariance(a,b){
 339   return polygamma(1,a)-polygamma(a+b);
 340 }
 341
 342 define betaskewness(a,b){
 343   return (2*(b-a)*sqrt(a+b+1))/( (a+b+1)*sqrt(a*b) );
 344 }
 345
 346 define betakurtosis(a,b){
 347   local num denom;
 348
 349   num = 6*( (a-b)^2*(a+b+1)-a*b*(a+b+2));
 350   denom = a*b*(a+b+2)*(a+b+3);
 351   return num/denom;
 352 }
 353
 354 define betaentropy(a,b){
 355   return lnbeta(a,b)-(a-1)*psi(a)-(b-1)*psi(b)+(a+b+1)*psi(a+b);
 356
 357 }
 358
 359 /*******************************************************************************
 360  *
 361  *
 362  *                Normal (Gaussian) distribution
 363  *
 364  *
 365  ******************************************************************************/
 366
 367
 368 define normalpdf(x,mu,sigma){
 369   return 1/(sqrt(2*pi()*sigma^2))*exp( ( (x-mu)^2 )/( 2*sigma^2 ) );
 370 }
 371
 372 define normalcdf(x,mu,sigma){
 373   return 1/2*(1+erf( ( x-mu )/( sqrt(2*sigma^2) )  )  );
 374 }
 375
 376 define probit(p){
 377   if(p<0 || p > 1) return newerror("probit: p out of domain 0<=p<=1");
 378   return sqrt(2)*ervinv(2*p-1);
 379 }
 380
 381 define normalcdfinv(p,mu,sigma){
 382   if(p<0 || p > 1) return newerror("normalcdfinv: p out of domain 0<=p<=1");
 383   return mu+ sigma*probit(p);
 384 }
 385
 386 define normalmean(mu,sigma){return mu;}
 387
 388 define normalmedian(mu,sigma){return mu;}
 389
 390 define normalmode(mu,sigma){return mu;}
 391
 392 define normalvariance(mu,sigma){return sigma^2;}
 393
 394 define normalskewness(mu,sigma){return 0;}
 395
 396 define normalkurtosis(mu,sigma){return 0;}
 397
 398 define normalentropy(mu,sigma){
 399   return 1/3*ln( 2*pi()*exp(1)*sigma^2 );
 400 }
 401
 402 /* moment generating f. */
 403 define normalmgf(mu,sigma,t){
 404   return exp(mu*t+1/2*sigma^2*t^2);
 405 }
 406
 407 /* characteristic f. */
 408 define normalcf(mu,sigma,t){
 409   return exp(mu*t-1/2*sigma^2*t^2);
 410 }
 411
 412
 413 /*******************************************************************************
 414  *
 415  *
 416  *                Chi-squared distribution
 417  *
 418  *
 419  ******************************************************************************/
 420
 421 define chisquaredpdf(x,k){
 422   if(!isint(k) || k<0) return newerror("chisquaredpdf: k not in N");
 423   if(im(x) || x<0) return newerror("chisquaredpdf: x not in +R");
 424   /* The gamma function does not check for half integers, do it here? */
 425   return 1/(2^(k/2)*gamma(k/2))*x^((k/2)-1)*exp(-x/2);
 426 }
 427
 428 define chisquaredpcdf(x,k){
 429   if(!isint(k) || k<0) return newerror("chisquaredcdf: k not in N");
 430   if(im(x) || x<0) return newerror("chisquaredcdf: x not in +R");
 431
 432   return 1/(gamma(k/2))*gammainc(k/2,x/2);
 433 }
 434
 435 define chisquaredmean(x,k){return k;}
 436
 437 define chisquaredmedian(x,k){
 438   /* TODO: implement a FAST inverse incomplete gamma-{q,p} function */
 439   return k*(1-2/(9*k))^3;
 440 }
 441
 442 define chisquaredmode(x,k){return max(k-2,0);}
 443 define chisquaredvariance(x,k){return 2*k;}
 444 define chisquaredskewness(x,k){return sqrt(8/k);}
 445 define chisquaredkurtosis(x,k){return 12/k;}
 446 define chisquaredentropy(x,k){
 447   return k/2+ln(2*gamma(k/2)) + (1-k/2)*psi(k/2);
 448 }
 449
 450 define chisquaredmfg(k,t){
 451   if(t>=1/2)return newerror("chisquaredmfg: t >= 1/2");
 452   return (1-2*t)^(k/2);
 453 }
 454
 455 define chisquaredcf(k,t){
 456   return (1-2*1i*t)^(k/2);
 457 }
 458
 459
 460 /*
 461  * restore internal function from resource debugging
 462  */
 463 config("resource_debug", resource_debug_level),;
 464 if (config("resource_debug") & 3) {
 465     print "gammaincoctave(z,a)";
 466     print "invbetainc(x,a,b)";
 467     print "betapdf(x,a,b)";
 468     print "betacdf(x,a,b)";
 469     print "betacdfinv(x,a,b)";
 470     print "betamedian(a,b)";
 471     print "betamode(a,b)";
 472     print "betavariance(a,b)";
 473     print "betalnvariance(a,b)";
 474     print "betaskewness(a,b)";
 475     print "betakurtosis(a,b)";
 476     print "betaentropy(a,b)";
 477     print "normalpdf(x,mu,sigma)";
 478     print "normalcdf(x,mu,sigma)";
 479     print "probit(p)";
 480     print "normalcdfinv(p,mu,sigma)";
 481     print "normalmean(mu,sigma)";
 482     print "normalmedian(mu,sigma)";
 483     print "normalmode(mu,sigma)";
 484     print "normalvariance(mu,sigma)";
 485     print "normalskewness(mu,sigma)";
 486     print "normalkurtosis(mu,sigma)";
 487     print "normalentropy(mu,sigma)";
 488     print "normalmgf(mu,sigma,t)";
 489     print "normalcf(mu,sigma,t)";
 490     print "chisquaredpdf(x,k)";
 491     print "chisquaredpcdf(x,k)";
 492     print "chisquaredmean(x,k)";
 493     print "chisquaredmedian(x,k)";
 494     print "chisquaredmode(x,k)";
 495     print "chisquaredvariance(x,k)";
 496     print "chisquaredskewness(x,k)";
 497     print "chisquaredkurtosis(x,k)";
 498     print "chisquaredentropy(x,k)";
 499     print "chisquaredmfg(k,t)";
 500     print "chisquaredcf(k,t)";
 501 }
 502