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