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