/*
Copyright � 1999 CERN - European Organization for Nuclear Research.
Permission to use, copy, modify, distribute and sell this software and its documentation for any purpose
is hereby granted without fee, provided that the above copyright notice appear in all copies and
that both that copyright notice and this permission notice appear in supporting documentation.
CERN makes no representations about the suitability of this software for any purpose.
It is provided "as is" without expressed or implied warranty.
*/
package sim.util.distribution;
import ec.util.MersenneTwisterFast;
/**
* Custom tailored numerical integration of certain probability distributions.
* <p>
* <b>Implementation:</b>
* <dt>
* Some code taken and adapted from the <A HREF="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">Java 2D Graph Package 2.4</A>,
* which in turn is a port from the <A HREF="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</A> Math Library (C).
* Most Cephes code (missing from the 2D Graph Package) directly ported.
*
* @author peter.gedeck@pharma.Novartis.com
* @author wolfgang.hoschek@cern.ch
* @version 0.91, 08-Dec-99
*/
public class Probability extends Constants {
private static final long serialVersionUID = 1;
/*************************************************
* COEFFICIENTS FOR METHOD normalInverse() *
*************************************************/
/* approximation for 0 <= |y - 0.5| <= 3/8 */
protected static final double P0[] = {
-5.99633501014107895267E1,
9.80010754185999661536E1,
-5.66762857469070293439E1,
1.39312609387279679503E1,
-1.23916583867381258016E0,
};
protected static final double Q0[] = {
/* 1.00000000000000000000E0,*/
1.95448858338141759834E0,
4.67627912898881538453E0,
8.63602421390890590575E1,
-2.25462687854119370527E2,
2.00260212380060660359E2,
-8.20372256168333339912E1,
1.59056225126211695515E1,
-1.18331621121330003142E0,
};
/* Approximation for interval z = sqrt(-2 log y ) between 2 and 8
* i.e., y between exp(-2) = .135 and exp(-32) = 1.27e-14.
*/
protected static final double P1[] = {
4.05544892305962419923E0,
3.15251094599893866154E1,
5.71628192246421288162E1,
4.40805073893200834700E1,
1.46849561928858024014E1,
2.18663306850790267539E0,
-1.40256079171354495875E-1,
-3.50424626827848203418E-2,
-8.57456785154685413611E-4,
};
protected static final double Q1[] = {
/* 1.00000000000000000000E0,*/
1.57799883256466749731E1,
4.53907635128879210584E1,
4.13172038254672030440E1,
1.50425385692907503408E1,
2.50464946208309415979E0,
-1.42182922854787788574E-1,
-3.80806407691578277194E-2,
-9.33259480895457427372E-4,
};
/* Approximation for interval z = sqrt(-2 log y ) between 8 and 64
* i.e., y between exp(-32) = 1.27e-14 and exp(-2048) = 3.67e-890.
*/
protected static final double P2[] = {
3.23774891776946035970E0,
6.91522889068984211695E0,
3.93881025292474443415E0,
1.33303460815807542389E0,
2.01485389549179081538E-1,
1.23716634817820021358E-2,
3.01581553508235416007E-4,
2.65806974686737550832E-6,
6.23974539184983293730E-9,
};
protected static final double Q2[] = {
/* 1.00000000000000000000E0,*/
6.02427039364742014255E0,
3.67983563856160859403E0,
1.37702099489081330271E0,
2.16236993594496635890E-1,
1.34204006088543189037E-2,
3.28014464682127739104E-4,
2.89247864745380683936E-6,
6.79019408009981274425E-9,
};
/**
* Makes this class non instantiable, but still let's others inherit from it.
*/
protected Probability() {}
/**
* Returns the area from zero to <tt>x</tt> under the beta density
* function.
* <pre>
* x
* - -
* | (a+b) | | a-1 b-1
* P(x) = ---------- | t (1-t) dt
* - - | |
* | (a) | (b) -
* 0
* </pre>
* This function is identical to the incomplete beta
* integral function <tt>incompleteBeta(a, b, x)</tt>.
*
* The complemented function is
*
* <tt>1 - P(1-x) = incompleteBeta( b, a, x )</tt>;
*
*/
static public double beta(double a, double b, double x ) {
return incompleteBeta( a, b, x );
}
/**
* Returns the area under the right hand tail (from <tt>x</tt> to
* infinity) of the beta density function.
*
* This function is identical to the incomplete beta
* integral function <tt>incompleteBeta(b, a, x)</tt>.
*/
static public double betaComplemented(double a, double b, double x ) {
return incompleteBeta( b, a, x );
}
/**
* Returns the sum of the terms <tt>0</tt> through <tt>k</tt> of the Binomial
* probability density.
* <pre>
* k
* -- ( n ) j n-j
* > ( ) p (1-p)
* -- ( j )
* j=0
* </pre>
* The terms are not summed directly; instead the incomplete
* beta integral is employed, according to the formula
* <p>
* <tt>y = binomial( k, n, p ) = incompleteBeta( n-k, k+1, 1-p )</tt>.
* <p>
* All arguments must be positive,
* @param k end term.
* @param n the number of trials.
* @param p the probability of success (must be in <tt>(0.0,1.0)</tt>).
*/
static public double binomial(int k, int n, double p) {
if( (p < 0.0) || (p > 1.0) ) throw new IllegalArgumentException();
if( (k < 0) || (n < k) ) throw new IllegalArgumentException();
if( k == n ) return( 1.0 );
if( k == 0 ) return Math.pow( 1.0-p, n-k );
return incompleteBeta( n-k, k+1, 1.0 - p );
}
/**
* Returns the sum of the terms <tt>k+1</tt> through <tt>n</tt> of the Binomial
* probability density.
* <pre>
* n
* -- ( n ) j n-j
* > ( ) p (1-p)
* -- ( j )
* j=k+1
* </pre>
* The terms are not summed directly; instead the incomplete
* beta integral is employed, according to the formula
* <p>
* <tt>y = binomialComplemented( k, n, p ) = incompleteBeta( k+1, n-k, p )</tt>.
* <p>
* All arguments must be positive,
* @param k end term.
* @param n the number of trials.
* @param p the probability of success (must be in <tt>(0.0,1.0)</tt>).
*/
static public double binomialComplemented(int k, int n, double p) {
if( (p < 0.0) || (p > 1.0) ) throw new IllegalArgumentException();
if( (k < 0) || (n < k) ) throw new IllegalArgumentException();
if( k == n ) return( 0.0 );
if( k == 0 ) return 1.0 - Math.pow( 1.0-p, n-k );
return incompleteBeta( k+1, n-k, p );
}
/**
* Returns the area under the left hand tail (from 0 to <tt>x</tt>)
* of the Chi square probability density function with
* <tt>v</tt> degrees of freedom.
* <pre>
* inf.
* -
* 1 | | v/2-1 -t/2
* P( x | v ) = ----------- | t e dt
* v/2 - | |
* 2 | (v/2) -
* x
* </pre>
* where <tt>x</tt> is the Chi-square variable.
* <p>
* The incomplete gamma integral is used, according to the
* formula
* <p>
* <tt>y = chiSquare( v, x ) = incompleteGamma( v/2.0, x/2.0 )</tt>.
* <p>
* The arguments must both be positive.
*
* @param v degrees of freedom.
* @param x integration end point.
*/
static public double chiSquare(double v, double x) throws ArithmeticException {
if( x < 0.0 || v < 1.0 ) return 0.0;
return incompleteGamma( v/2.0, x/2.0 );
}
/**
* Returns the area under the right hand tail (from <tt>x</tt> to
* infinity) of the Chi square probability density function
* with <tt>v</tt> degrees of freedom.
* <pre>
* inf.
* -
* 1 | | v/2-1 -t/2
* P( x | v ) = ----------- | t e dt
* v/2 - | |
* 2 | (v/2) -
* x
* </pre>
* where <tt>x</tt> is the Chi-square variable.
*
* The incomplete gamma integral is used, according to the
* formula
*
* <tt>y = chiSquareComplemented( v, x ) = incompleteGammaComplement( v/2.0, x/2.0 )</tt>.
*
*
* The arguments must both be positive.
*
* @param v degrees of freedom.
*/
static public double chiSquareComplemented(double v, double x) throws ArithmeticException {
if( x < 0.0 || v < 1.0 ) return 0.0;
return incompleteGammaComplement( v/2.0, x/2.0 );
}
/**
* Returns the error function of the normal distribution; formerly named <tt>erf</tt>.
* The integral is
* <pre>
* x
* -
* 2 | | 2
* erf(x) = -------- | exp( - t ) dt.
* sqrt(pi) | |
* -
* 0
* </pre>
* <b>Implementation:</b>
* For <tt>0 <= |x| < 1, erf(x) = x * P4(x**2)/Q5(x**2)</tt>; otherwise
* <tt>erf(x) = 1 - erfc(x)</tt>.
* <p>
* Code adapted from the <A HREF="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">Java 2D Graph Package 2.4</A>,
* which in turn is a port from the <A HREF="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</A> Math Library (C).
*
* @param a the argument to the function.
*/
static public double errorFunction(double x) throws ArithmeticException {
double y, z;
final double T[] = {
9.60497373987051638749E0,
9.00260197203842689217E1,
2.23200534594684319226E3,
7.00332514112805075473E3,
5.55923013010394962768E4
};
final double U[] = {
//1.00000000000000000000E0,
3.35617141647503099647E1,
5.21357949780152679795E2,
4.59432382970980127987E3,
2.26290000613890934246E4,
4.92673942608635921086E4
};
if( Math.abs(x) > 1.0 ) return( 1.0 - errorFunctionComplemented(x) );
z = x * x;
y = x * Polynomial.polevl( z, T, 4 ) / Polynomial.p1evl( z, U, 5 );
return y;
}
/**
* Returns the complementary Error function of the normal distribution; formerly named <tt>erfc</tt>.
* <pre>
* 1 - erf(x) =
*
* inf.
* -
* 2 | | 2
* erfc(x) = -------- | exp( - t ) dt
* sqrt(pi) | |
* -
* x
* </pre>
* <b>Implementation:</b>
* For small x, <tt>erfc(x) = 1 - erf(x)</tt>; otherwise rational
* approximations are computed.
* <p>
* Code adapted from the <A HREF="http://www.sci.usq.edu.au/staff/leighb/graph/Top.html">Java 2D Graph Package 2.4</A>,
* which in turn is a port from the <A HREF="http://people.ne.mediaone.net/moshier/index.html#Cephes">Cephes 2.2</A> Math Library (C).
*
* @param a the argument to the function.
*/
static public double errorFunctionComplemented(double a) throws ArithmeticException {
double x,y,z,p,q;
double P[] = {
2.46196981473530512524E-10,
5.64189564831068821977E-1,
7.46321056442269912687E0,
4.86371970985681366614E1,
1.96520832956077098242E2,
5.26445194995477358631E2,
9.34528527171957607540E2,
1.02755188689515710272E3,
5.57535335369399327526E2
};
double Q[] = {
//1.0
1.32281951154744992508E1,
8.67072140885989742329E1,
3.54937778887819891062E2,
9.75708501743205489753E2,
1.82390916687909736289E3,
2.24633760818710981792E3,
1.65666309194161350182E3,
5.57535340817727675546E2
};
double R[] = {
5.64189583547755073984E-1,
1.27536670759978104416E0,
5.01905042251180477414E0,
6.16021097993053585195E0,
7.40974269950448939160E0,
2.97886665372100240670E0
};
double S[] = {
//1.00000000000000000000E0,
2.26052863220117276590E0,
9.39603524938001434673E0,
1.20489539808096656605E1,
1.70814450747565897222E1,
9.60896809063285878198E0,
3.36907645100081516050E0
};
if( a < 0.0 ) x = -a;
else x = a;
if( x < 1.0 ) return 1.0 - errorFunction(a);
z = -a * a;
if( z < -MAXLOG ) {
if( a < 0 ) return( 2.0 );
else return( 0.0 );
}
z = Math.exp(z);
if( x < 8.0 ) {
p = Polynomial.polevl( x, P, 8 );
q = Polynomial.p1evl( x, Q, 8 );
} else {
p = Polynomial.polevl( x, R, 5 );
q = Polynomial.p1evl( x, S, 6 );
}
y = (z * p)/q;
if( a < 0 ) y = 2.0 - y;
if( y == 0.0 ) {
if( a < 0 ) return 2.0;
else return( 0.0 );
}
return y;
}
/**
* Returns the integral from zero to <tt>x</tt> of the gamma probability
* density function.
* <pre>
* x
* b -
* a | | b-1 -at
* y = ----- | t e dt
* - | |
* | (b) -
* 0
* </pre>
* The incomplete gamma integral is used, according to the
* relation
*
* <tt>y = incompleteGamma( b, a*x )</tt>.
*
* @param a the paramater a (alpha) of the gamma distribution.
* @param b the paramater b (beta, lambda) of the gamma distribution.
* @param x integration end point.
*/
static public double gamma(double a, double b, double x ) {
if( x < 0.0 ) return 0.0;
return incompleteGamma(b, a*x);
}
/**
* Returns the integral from <tt>x</tt> to infinity of the gamma
* probability density function:
* <pre>
* inf.
* b -
* a | | b-1 -at
* y = ----- | t e dt
* - | |
* | (b) -
* x
* </pre>
* The incomplete gamma integral is used, according to the
* relation
* <p>
* y = incompleteGammaComplement( b, a*x ).
*
* @param a the paramater a (alpha) of the gamma distribution.
* @param b the paramater b (beta, lambda) of the gamma distribution.
* @param x integration end point.
*/
static public double gammaComplemented(double a, double b, double x ) {
if( x < 0.0 ) return 0.0;
return incompleteGammaComplement(b, a*x);
}
/**
* Returns the sum of the terms <tt>0</tt> through <tt>k</tt> of the Negative Binomial Distribution.
* <pre>
* k
* -- ( n+j-1 ) n j
* > ( ) p (1-p)
* -- ( j )
* j=0
* </pre>
* In a sequence of Bernoulli trials, this is the probability
* that <tt>k</tt> or fewer failures precede the <tt>n</tt>-th success.
* <p>
* The terms are not computed individually; instead the incomplete
* beta integral is employed, according to the formula
* <p>
* <tt>y = negativeBinomial( k, n, p ) = incompleteBeta( n, k+1, p )</tt>.
*
* All arguments must be positive,
* @param k end term.
* @param n the number of trials.
* @param p the probability of success (must be in <tt>(0.0,1.0)</tt>).
*/
static public double negativeBinomial(int k, int n, double p) {
if( (p < 0.0) || (p > 1.0) ) throw new IllegalArgumentException();
if(k < 0) return 0.0;
return incompleteBeta( n, k+1, p );
}
/**
* Returns the sum of the terms <tt>k+1</tt> to infinity of the Negative
* Binomial distribution.
* <pre>
* inf
* -- ( n+j-1 ) n j
* > ( ) p (1-p)
* -- ( j )
* j=k+1
* </pre>
* The terms are not computed individually; instead the incomplete
* beta integral is employed, according to the formula
* <p>
* y = negativeBinomialComplemented( k, n, p ) = incompleteBeta( k+1, n, 1-p ).
*
* All arguments must be positive,
* @param k end term.
* @param n the number of trials.
* @param p the probability of success (must be in <tt>(0.0,1.0)</tt>).
*/
static public double negativeBinomialComplemented(int k, int n, double p) {
if( (p < 0.0) || (p > 1.0) ) throw new IllegalArgumentException();
if(k < 0) return 0.0;
return incompleteBeta( k+1, n, 1.0-p );
}
/**
* Returns the area under the Normal (Gaussian) probability density
* function, integrated from minus infinity to <tt>x</tt> (assumes mean is zero, variance is one).
* <pre>
* x
* -
* 1 | | 2
* normal(x) = --------- | exp( - t /2 ) dt
* sqrt(2pi) | |
* -
* -inf.
*
* = ( 1 + erf(z) ) / 2
* = erfc(z) / 2
* </pre>
* where <tt>z = x/sqrt(2)</tt>.
* Computation is via the functions <tt>errorFunction</tt> and <tt>errorFunctionComplement</tt>.
*/
static public double normal( double a) throws ArithmeticException {
double x, y, z;
x = a * SQRTH;
z = Math.abs(x);
if( z < SQRTH ) y = 0.5 + 0.5 * errorFunction(x);
else {
y = 0.5 * errorFunctionComplemented(z);
if( x > 0 ) y = 1.0 - y;
}
return y;
}
/**
* Returns the area under the Normal (Gaussian) probability density
* function, integrated from minus infinity to <tt>x</tt>.
* <pre>
* x
* -
* 1 | | 2
* normal(x) = --------- | exp( - (t-mean) / 2v ) dt
* sqrt(2pi*v)| |
* -
* -inf.
*
* </pre>
* where <tt>v = variance</tt>.
* Computation is via the functions <tt>errorFunction</tt>.
*
* @param mean the mean of the normal distribution.
* @param variance the variance of the normal distribution.
* @param x the integration limit.
*/
static public double normal(double mean, double variance, double x) throws ArithmeticException {
if (x>0)
return 0.5 + 0.5*errorFunction((x-mean)/Math.sqrt(2.0*variance));
else
return 0.5 - 0.5*errorFunction((-(x-mean))/Math.sqrt(2.0*variance));
}
/**
* Returns the value, <tt>x</tt>, for which the area under the
* Normal (Gaussian) probability density function (integrated from
* minus infinity to <tt>x</tt>) is equal to the argument <tt>y</tt> (assumes mean is zero, variance is one); formerly named <tt>ndtri</tt>.
* <p>
* For small arguments <tt>0 < y < exp(-2)</tt>, the program computes
* <tt>z = sqrt( -2.0 * log(y) )</tt>; then the approximation is
* <tt>x = z - log(z)/z - (1/z) P(1/z) / Q(1/z)</tt>.
* There are two rational functions P/Q, one for <tt>0 < y < exp(-32)</tt>
* and the other for <tt>y</tt> up to <tt>exp(-2)</tt>.
* For larger arguments,
* <tt>w = y - 0.5</tt>, and <tt>x/sqrt(2pi) = w + w**3 R(w**2)/S(w**2))</tt>.
*
*/
static public double normalInverse( double y0) throws ArithmeticException {
double x, y, z, y2, x0, x1;
int code;
final double s2pi = Math.sqrt(2.0*Math.PI);
if( y0 <= 0.0 ) throw new IllegalArgumentException();
if( y0 >= 1.0 ) throw new IllegalArgumentException();
code = 1;
y = y0;
if( y > (1.0 - 0.13533528323661269189) ) { /* 0.135... = exp(-2) */
y = 1.0 - y;
code = 0;
}
if( y > 0.13533528323661269189 ) {
y = y - 0.5;
y2 = y * y;
x = y + y * (y2 * Polynomial.polevl( y2, P0, 4)/Polynomial.p1evl( y2, Q0, 8 ));
x = x * s2pi;
return(x);
}
x = Math.sqrt( -2.0 * Math.log(y) );
x0 = x - Math.log(x)/x;
z = 1.0/x;
if( x < 8.0 ) /* y > exp(-32) = 1.2664165549e-14 */
x1 = z * Polynomial.polevl( z, P1, 8 )/Polynomial.p1evl( z, Q1, 8 );
else
x1 = z * Polynomial.polevl( z, P2, 8 )/Polynomial.p1evl( z, Q2, 8 );
x = x0 - x1;
if( code != 0 )
x = -x;
return( x );
}
/**
* Returns the sum of the first <tt>k</tt> terms of the Poisson distribution.
* <pre>
* k j
* -- -m m
* > e --
* -- j!
* j=0
* </pre>
* The terms are not summed directly; instead the incomplete
* gamma integral is employed, according to the relation
* <p>
* <tt>y = poisson( k, m ) = incompleteGammaComplement( k+1, m )</tt>.
*
* The arguments must both be positive.
*
* @param k number of terms.
* @param mean the mean of the poisson distribution.
*/
static public double poisson(int k, double mean) throws ArithmeticException {
if( mean < 0 ) throw new IllegalArgumentException();
if( k < 0 ) return 0.0;
return incompleteGammaComplement((double)(k+1) ,mean);
}
/**
* Returns the sum of the terms <tt>k+1</tt> to <tt>Infinity</tt> of the Poisson distribution.
* <pre>
* inf. j
* -- -m m
* > e --
* -- j!
* j=k+1
* </pre>
* The terms are not summed directly; instead the incomplete
* gamma integral is employed, according to the formula
* <p>
* <tt>y = poissonComplemented( k, m ) = incompleteGamma( k+1, m )</tt>.
*
* The arguments must both be positive.
*
* @param k start term.
* @param mean the mean of the poisson distribution.
*/
static public double poissonComplemented(int k, double mean) throws ArithmeticException {
if( mean < 0 ) throw new IllegalArgumentException();
if( k < -1 ) return 0.0;
return incompleteGamma((double)(k+1),mean);
}
/**
* Returns the integral from minus infinity to <tt>t</tt> of the Student-t
* distribution with <tt>k > 0</tt> degrees of freedom.
* <pre>
* t
* -
* | |
* - | 2 -(k+1)/2
* | ( (k+1)/2 ) | ( x )
* ---------------------- | ( 1 + --- ) dx
* - | ( k )
* sqrt( k pi ) | ( k/2 ) |
* | |
* -
* -inf.
* </pre>
* Relation to incomplete beta integral:
* <p>
* <tt>1 - studentT(k,t) = 0.5 * incompleteBeta( k/2, 1/2, z )</tt>
* where <tt>z = k/(k + t**2)</tt>.
* <p>
* Since the function is symmetric about <tt>t=0</tt>, the area under the
* right tail of the density is found by calling the function
* with <tt>-t</tt> instead of <tt>t</tt>.
*
* @param k degrees of freedom.
* @param t integration end point.
*/
static public double studentT(double k, double t) throws ArithmeticException {
if( k <= 0 ) throw new IllegalArgumentException();
if( t == 0 ) return( 0.5 );
double cdf = 0.5 * incompleteBeta( 0.5*k, 0.5, k / (k + t * t) );
if (t >= 0) cdf = 1.0 - cdf; // fixes bug reported by stefan.bentink@molgen.mpg.de
return cdf;
}
/**
* Returns the value, <tt>t</tt>, for which the area under the
* Student-t probability density function (integrated from
* minus infinity to <tt>t</tt>) is equal to <tt>1-alpha/2</tt>.
* The value returned corresponds to usual Student t-distribution lookup
* table for <tt>t<sub>alpha[size]</sub></tt>.
* <p>
* The function uses the studentT function to determine the return
* value iteratively.
*
* @param alpha probability
* @param size size of data set
*/
public static double studentTInverse(double alpha, int size) {
double cumProb = 1-alpha/2; // Cumulative probability
double f1,f2,f3;
double x1,x2,x3;
double g,s12;
cumProb = 1-alpha/2; // Cumulative probability
x1 = normalInverse(cumProb);
// Return inverse of normal for large size
if (size > 200) {
return x1;
}
// Find a pair of x1,x2 that braket zero
f1 = studentT(size,x1)-cumProb;
x2 = x1; f2 = f1;
do {
if (f1>0) {
x2 = x2/2;
} else {
x2 = x2+x1;
}
f2 = studentT(size,x2)-cumProb;
} while (f1*f2>0);
// Find better approximation
// Pegasus-method
do {
// Calculate slope of secant and t value for which it is 0.
s12 = (f2-f1)/(x2-x1);
x3 = x2 - f2/s12;
// Calculate function value at x3
f3 = studentT(size,x3)-cumProb;
if (Math.abs(f3)<1e-8) { // This criteria needs to be very tight!
// We found a perfect value -> return
return x3;
}
if (f3*f2<0) {
x1=x2; f1=f2;
x2=x3; f2=f3;
} else {
g = f2/(f2+f3);
f1=g*f1;
x2=x3; f2=f3;
}
} while(Math.abs(x2-x1)>0.001);
if (Math.abs(f2)<=Math.abs(f1)) {
return x2;
} else {
return x1;
}
}
/////// MERGED IN BY SEAN FROM stat.Gamma
/**
* Returns the beta function of the arguments.
* <pre>
* - -
* | (a) | (b)
* beta( a, b ) = -----------.
* -
* | (a+b)
* </pre>
*/
static public double beta(double a, double b) throws ArithmeticException {
double y;
y = a + b;
y = gamma(y);
if( y == 0.0 ) return 1.0;
if( a > b ) {
y = gamma(a)/y;
y *= gamma(b);
}
else {
y = gamma(b)/y;
y *= gamma(a);
}
return(y);
}
/**
* Returns the Gamma function of the argument.
*/
static public double gamma(double x) throws ArithmeticException {
double P[] = {
1.60119522476751861407E-4,
1.19135147006586384913E-3,
1.04213797561761569935E-2,
4.76367800457137231464E-2,
2.07448227648435975150E-1,
4.94214826801497100753E-1,
9.99999999999999996796E-1
};
double Q[] = {
-2.31581873324120129819E-5,
5.39605580493303397842E-4,
-4.45641913851797240494E-3,
1.18139785222060435552E-2,
3.58236398605498653373E-2,
-2.34591795718243348568E-1,
7.14304917030273074085E-2,
1.00000000000000000320E0
};
//double MAXGAM = 171.624376956302725;
//double LOGPI = 1.14472988584940017414;
double p, z;
int i;
double q = Math.abs(x);
if( q > 33.0 ) {
if( x < 0.0 ) {
p = Math.floor(q);
if( p == q ) throw new ArithmeticException("gamma: overflow");
//i = (int)p;
z = q - p;
if( z > 0.5 ) {
p += 1.0;
z = q - p;
}
z = q * Math.sin( Math.PI * z );
if( z == 0.0 ) throw new ArithmeticException("gamma: overflow");
z = Math.abs(z);
z = Math.PI/(z * stirlingFormula(q) );
return -z;
} else {
return stirlingFormula(x);
}
}
z = 1.0;
while( x >= 3.0 ) {
x -= 1.0;
z *= x;
}
while( x < 0.0 ) {
if( x == 0.0 ) {
throw new ArithmeticException("gamma: singular");
} else
if( x > -1.E-9 ) {
return( z/((1.0 + 0.5772156649015329 * x) * x) );
}
z /= x;
x += 1.0;
}
while( x < 2.0 ) {
if( x == 0.0 ) {
throw new ArithmeticException("gamma: singular");
} else
if( x < 1.e-9 ) {
return( z/((1.0 + 0.5772156649015329 * x) * x) );
}
z /= x;
x += 1.0;
}
if( (x == 2.0) || (x == 3.0) ) return z;
x -= 2.0;
p = Polynomial.polevl( x, P, 6 );
q = Polynomial.polevl( x, Q, 7 );
return z * p / q;
}
/**
* Returns the Incomplete Beta Function evaluated from zero to <tt>xx</tt>; formerly named <tt>ibeta</tt>.
*
* @param aa the alpha parameter of the beta distribution.
* @param bb the beta parameter of the beta distribution.
* @param xx the integration end point.
*/
public static double incompleteBeta( double aa, double bb, double xx ) throws ArithmeticException {
double a, b, t, x, xc, w, y;
boolean flag;
if( aa <= 0.0 || bb <= 0.0 ) throw new
ArithmeticException("ibeta: Domain error!");
if( (xx <= 0.0) || ( xx >= 1.0) ) {
if( xx == 0.0 ) return 0.0;
if( xx == 1.0 ) return 1.0;
throw new ArithmeticException("ibeta: Domain error!");
}
flag = false;
if( (bb * xx) <= 1.0 && xx <= 0.95) {
t = powerSeries(aa, bb, xx);
return t;
}
w = 1.0 - xx;
/* Reverse a and b if x is greater than the mean. */
if( xx > (aa/(aa+bb)) ) {
flag = true;
a = bb;
b = aa;
xc = xx;
x = w;
} else {
a = aa;
b = bb;
xc = w;
x = xx;
}
if( flag && (b * x) <= 1.0 && x <= 0.95) {
t = powerSeries(a, b, x);
if( t <= MACHEP ) t = 1.0 - MACHEP;
else t = 1.0 - t;
return t;
}
/* Choose expansion for better convergence. */
y = x * (a+b-2.0) - (a-1.0);
if( y < 0.0 )
w = incompleteBetaFraction1( a, b, x );
else
w = incompleteBetaFraction2( a, b, x ) / xc;
/* Multiply w by the factor
a b _ _ _
x (1-x) | (a+b) / ( a | (a) | (b) ) . */
y = a * Math.log(x);
t = b * Math.log(xc);
if( (a+b) < MAXGAM && Math.abs(y) < MAXLOG && Math.abs(t) < MAXLOG ) {
t = Math.pow(xc,b);
t *= Math.pow(x,a);
t /= a;
t *= w;
t *= gamma(a+b) / (gamma(a) * gamma(b));
if( flag ) {
if( t <= MACHEP ) t = 1.0 - MACHEP;
else t = 1.0 - t;
}
return t;
}
/* Resort to logarithms. */
y += t + logGamma(a+b) - logGamma(a) - logGamma(b);
y += Math.log(w/a);
if( y < MINLOG )
t = 0.0;
else
t = Math.exp(y);
if( flag ) {
if( t <= MACHEP ) t = 1.0 - MACHEP;
else t = 1.0 - t;
}
return t;
}
/**
* Continued fraction expansion #1 for incomplete beta integral; formerly named <tt>incbcf</tt>.
*/
static double incompleteBetaFraction1( double a, double b, double x ) throws ArithmeticException {
double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
double k1, k2, k3, k4, k5, k6, k7, k8;
double r, t, ans, thresh;
int n;
k1 = a;
k2 = a + b;
k3 = a;
k4 = a + 1.0;
k5 = 1.0;
k6 = b - 1.0;
k7 = k4;
k8 = a + 2.0;
pkm2 = 0.0;
qkm2 = 1.0;
pkm1 = 1.0;
qkm1 = 1.0;
ans = 1.0;
r = 1.0;
n = 0;
thresh = 3.0 * MACHEP;
do {
xk = -( x * k1 * k2 )/( k3 * k4 );
pk = pkm1 + pkm2 * xk;
qk = qkm1 + qkm2 * xk;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
xk = ( x * k5 * k6 )/( k7 * k8 );
pk = pkm1 + pkm2 * xk;
qk = qkm1 + qkm2 * xk;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
if( qk != 0 ) r = pk/qk;
if( r != 0 ) {
t = Math.abs( (ans - r)/r );
ans = r;
} else
t = 1.0;
if( t < thresh ) return ans;
k1 += 1.0;
k2 += 1.0;
k3 += 2.0;
k4 += 2.0;
k5 += 1.0;
k6 -= 1.0;
k7 += 2.0;
k8 += 2.0;
if( (Math.abs(qk) + Math.abs(pk)) > big ) {
pkm2 *= biginv;
pkm1 *= biginv;
qkm2 *= biginv;
qkm1 *= biginv;
}
if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {
pkm2 *= big;
pkm1 *= big;
qkm2 *= big;
qkm1 *= big;
}
} while( ++n < 300 );
return ans;
}
/**
* Continued fraction expansion #2 for incomplete beta integral; formerly named <tt>incbd</tt>.
*/
static double incompleteBetaFraction2( double a, double b, double x ) throws ArithmeticException {
double xk, pk, pkm1, pkm2, qk, qkm1, qkm2;
double k1, k2, k3, k4, k5, k6, k7, k8;
double r, t, ans, z, thresh;
int n;
k1 = a;
k2 = b - 1.0;
k3 = a;
k4 = a + 1.0;
k5 = 1.0;
k6 = a + b;
k7 = a + 1.0;;
k8 = a + 2.0;
pkm2 = 0.0;
qkm2 = 1.0;
pkm1 = 1.0;
qkm1 = 1.0;
z = x / (1.0-x);
ans = 1.0;
r = 1.0;
n = 0;
thresh = 3.0 * MACHEP;
do {
xk = -( z * k1 * k2 )/( k3 * k4 );
pk = pkm1 + pkm2 * xk;
qk = qkm1 + qkm2 * xk;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
xk = ( z * k5 * k6 )/( k7 * k8 );
pk = pkm1 + pkm2 * xk;
qk = qkm1 + qkm2 * xk;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
if( qk != 0 ) r = pk/qk;
if( r != 0 ) {
t = Math.abs( (ans - r)/r );
ans = r;
} else
t = 1.0;
if( t < thresh ) return ans;
k1 += 1.0;
k2 -= 1.0;
k3 += 2.0;
k4 += 2.0;
k5 += 1.0;
k6 += 1.0;
k7 += 2.0;
k8 += 2.0;
if( (Math.abs(qk) + Math.abs(pk)) > big ) {
pkm2 *= biginv;
pkm1 *= biginv;
qkm2 *= biginv;
qkm1 *= biginv;
}
if( (Math.abs(qk) < biginv) || (Math.abs(pk) < biginv) ) {
pkm2 *= big;
pkm1 *= big;
qkm2 *= big;
qkm1 *= big;
}
} while( ++n < 300 );
return ans;
}
/**
* Returns the Incomplete Gamma function; formerly named <tt>igamma</tt>.
* @param a the parameter of the gamma distribution.
* @param x the integration end point.
*/
static public double incompleteGamma(double a, double x)
throws ArithmeticException {
double ans, ax, c, r;
if( x <= 0 || a <= 0 ) return 0.0;
if( x > 1.0 && x > a ) return 1.0 - incompleteGammaComplement(a,x);
/* Compute x**a * exp(-x) / gamma(a) */
ax = a * Math.log(x) - x - logGamma(a);
if( ax < -MAXLOG ) return( 0.0 );
ax = Math.exp(ax);
/* power series */
r = a;
c = 1.0;
ans = 1.0;
do {
r += 1.0;
c *= x/r;
ans += c;
}
while( c/ans > MACHEP );
return( ans * ax/a );
}
/**
* Returns the Complemented Incomplete Gamma function; formerly named <tt>igamc</tt>.
* @param a the parameter of the gamma distribution.
* @param x the integration start point.
*/
static public double incompleteGammaComplement( double a, double x ) throws ArithmeticException {
double ans, ax, c, yc, r, t, y, z;
double pk, pkm1, pkm2, qk, qkm1, qkm2;
if( x <= 0 || a <= 0 ) return 1.0;
if( x < 1.0 || x < a ) return 1.0 - incompleteGamma(a,x);
ax = a * Math.log(x) - x - logGamma(a);
if( ax < -MAXLOG ) return 0.0;
ax = Math.exp(ax);
/* continued fraction */
y = 1.0 - a;
z = x + y + 1.0;
c = 0.0;
pkm2 = 1.0;
qkm2 = x;
pkm1 = x + 1.0;
qkm1 = z * x;
ans = pkm1/qkm1;
do {
c += 1.0;
y += 1.0;
z += 2.0;
yc = y * c;
pk = pkm1 * z - pkm2 * yc;
qk = qkm1 * z - qkm2 * yc;
if( qk != 0 ) {
r = pk/qk;
t = Math.abs( (ans - r)/r );
ans = r;
} else
t = 1.0;
pkm2 = pkm1;
pkm1 = pk;
qkm2 = qkm1;
qkm1 = qk;
if( Math.abs(pk) > big ) {
pkm2 *= biginv;
pkm1 *= biginv;
qkm2 *= biginv;
qkm1 *= biginv;
}
} while( t > MACHEP );
return ans * ax;
}
/**
* Returns the natural logarithm of the gamma function; formerly named <tt>lgamma</tt>.
*/
public static double logGamma(double x) throws ArithmeticException {
double p, q, w, z;
double A[] = {
8.11614167470508450300E-4,
-5.95061904284301438324E-4,
7.93650340457716943945E-4,
-2.77777777730099687205E-3,
8.33333333333331927722E-2
};
double B[] = {
-1.37825152569120859100E3,
-3.88016315134637840924E4,
-3.31612992738871184744E5,
-1.16237097492762307383E6,
-1.72173700820839662146E6,
-8.53555664245765465627E5
};
double C[] = {
/* 1.00000000000000000000E0, */
-3.51815701436523470549E2,
-1.70642106651881159223E4,
-2.20528590553854454839E5,
-1.13933444367982507207E6,
-2.53252307177582951285E6,
-2.01889141433532773231E6
};
if( x < -34.0 ) {
q = -x;
w = logGamma(q);
p = Math.floor(q);
if( p == q ) throw new ArithmeticException("lgam: Overflow");
z = q - p;
if( z > 0.5 ) {
p += 1.0;
z = p - q;
}
z = q * Math.sin( Math.PI * z );
if( z == 0.0 ) throw new
ArithmeticException("lgamma: Overflow");
z = LOGPI - Math.log( z ) - w;
return z;
}
if( x < 13.0 ) {
z = 1.0;
while( x >= 3.0 ) {
x -= 1.0;
z *= x;
}
while( x < 2.0 ) {
if( x == 0.0 ) throw new
ArithmeticException("lgamma: Overflow");
z /= x;
x += 1.0;
}
if( z < 0.0 ) z = -z;
if( x == 2.0 ) return Math.log(z);
x -= 2.0;
p = x * Polynomial.polevl( x, B, 5 ) / Polynomial.p1evl( x, C, 6);
return( Math.log(z) + p );
}
if( x > 2.556348e305 ) throw new
ArithmeticException("lgamma: Overflow");
q = ( x - 0.5 ) * Math.log(x) - x + 0.91893853320467274178;
//if( x > 1.0e8 ) return( q );
if( x > 1.0e8 ) return( q );
p = 1.0/(x*x);
if( x >= 1000.0 )
q += (( 7.9365079365079365079365e-4 * p
- 2.7777777777777777777778e-3) *p
+ 0.0833333333333333333333) / x;
else
q += Polynomial.polevl( p, A, 4 ) / x;
return q;
}
/**
* Power series for incomplete beta integral; formerly named <tt>pseries</tt>.
* Use when b*x is small and x not too close to 1.
*/
static double powerSeries( double a, double b, double x ) throws ArithmeticException {
double s, t, u, v, n, t1, z, ai;
ai = 1.0 / a;
u = (1.0 - b) * x;
v = u / (a + 1.0);
t1 = v;
t = u;
n = 2.0;
s = 0.0;
z = MACHEP * ai;
while( Math.abs(v) > z ) {
u = (n - b) * x / n;
t *= u;
v = t / (a + n);
s += v;
n += 1.0;
}
s += t1;
s += ai;
u = a * Math.log(x);
if( (a+b) < MAXGAM && Math.abs(u) < MAXLOG ) {
t = gamma(a+b)/(gamma(a)*gamma(b));
s = s * t * Math.pow(x,a);
} else {
t = logGamma(a+b) - logGamma(a) - logGamma(b) + u + Math.log(s);
if( t < MINLOG ) s = 0.0;
else s = Math.exp(t);
}
return s;
}
/**
* Returns the Gamma function computed by Stirling's formula; formerly named <tt>stirf</tt>.
* The polynomial STIR is valid for 33 <= x <= 172.
*/
static double stirlingFormula(double x) throws ArithmeticException {
double STIR[] = {
7.87311395793093628397E-4,
-2.29549961613378126380E-4,
-2.68132617805781232825E-3,
3.47222221605458667310E-3,
8.33333333333482257126E-2,
};
double MAXSTIR = 143.01608;
double w = 1.0/x;
double y = Math.exp(x);
w = 1.0 + w * Polynomial.polevl( w, STIR, 4 );
if( x > MAXSTIR ) {
/* Avoid overflow in Math.pow() */
double v = Math.pow( x, 0.5 * x - 0.25 );
y = v * (v / y);
} else {
y = Math.pow( x, x - 0.5 ) / y;
}
y = SQTPI * y * w;
return y;
}
}