/**
* Class probLib deals with evaluating PDFs at specific data points and
* parameter values. Additional distributions can be added.
*
* @author Eric Ward
* @version August 31, 2006
*/
public class probabilityLibrary
{
// these are constants
final double PI = 3.14159265358979;
final double ROOT2PI = 2.506628274631;
final double LOGRT2PI = 0.91893853320467;
final double LN2PI = 1.837877066409;
final double LNPI = 1.1447298858494;
final double LN2 = 0.693147180559945;
// these variables are used for the gammln function
private double f;
// these variables are used for the gammln function
private double x, y, tmp, ser;
private int j;
final double cof[]={76.18009172947146,-86.50532032941677,24.01409824083091,-1.231739572450155,0.1208650973866179e-2,-0.5395239384953e-5};
double ret;
double z, z2;
/**
* Constructor for objects of class probLib
*/
public probabilityLibrary()
{
}
/**
* Method probLib returns the log of the probability density function (PDF)
* for a wide range of distributions.
*
* @param x, the data
* @param distID, the integer representation of the distribution
* @param distType, whether the distribution is discrete (0) or not (1)
* @param P0, the first parameter of the distribution
* @param P1, the second parameter of the distribution
* @param P2, the third parameter of the distribution
* @return double rnd
*/
public double probLib(double x, int distID, int distType, double P0, double P1, double P2)
{
//x is the value to evaluate the function at
//distID is the ID of the distribution
//distType: 0 for continuous, 1 for discrete
//P is the parameters of the distribution, e.g. (P0 = mu, P1 = sigma) for normal
ret = 0.0;
if(distType == 0) {
switch(distID) {
case 1:
// normal distribution (validated)
// PARS[0] = mu, PARS[1] = sigma
// generate a normal(0,1) random variable
ret = -0.5*LN2PI -Math.log(P1) -0.5*(x - P0)*(x - P0)/(P1 * P1);
break;
case 2:
// continuous uniform (validated)
// PARS[0] = low, PARS[1] = high
ret = -Math.log(P1 - P0);
break;
case 3:
// log-uniform: p.d.f = 1.0 / (x * (b - a))
// PARS[0] = low, PARS[1] = high
ret = -Math.log(Math.exp(P1) - Math.exp(P0)) - Math.log(x);
break;
case 4:
// gamma distribution (validated)
// generate a gamma(a,b) random variable
ret = -gammln(P0) -P0*Math.log(P1) + (P0-1)*Math.log(x) -x/P1;
break;
case 5:
// centralized gamma distribution (validated)
// generate a gamma(a,b) random variable, de-mean it, and add location param
// subtract location parameter, add the mean
// EW commented out the following line to change to standard C-syntax (++) 03/07/06
//x = x - P2 + P0*P1;
x += (-P2 + P0*P1);
ret = -gammln(P0) -P0*Math.log(P1) + (P0-1)*Math.log(x) -x/P1;
break;
case 6:
// lognormal distribution (validated)
z = Math.log(x);
ret = -0.5*LN2PI -Math.log(P1) -0.5*(z - P0)*(z - P0)/(P1 * P1) - z;
break;
case 7:
// exponential (validated)
ret = -Math.log(P0) - x/P0;
break;
case 8:
// inverse gamma (validated)
ret = -gammln(P0) + P0*Math.log(P1) - (P0-1)*Math.log(x) - P1/x;
break;
case 9:
// logistic (validated)
z = (x - P0)/P1;
ret = -Math.log(P1) - z - 2.0*Math.log(1 + Math.exp(-z));
break;
case 10:
// beta (validated)
ret = -Math.log(betaFunction((int)P0,(int)P1)) + (P0-1)*Math.log(x) + (P1-1)*Math.log(1-x);
break;
case 11:
// f distribution (validated)
z = P0/P1;
z2 = 0.5*(P0+P1);
ret = gammln(z2) - gammln(0.5*P0) - gammln(0.5*P1) + (0.5*P0)*Math.log(z) +
(0.5*(P0-2.0))*Math.log(x) - (z2)*Math.log(1 + x*z);
break;
case 12:
// chi-square
z = 0.5*P0;
ret = -gammln(z) -(z)*LN2 + (z-1)*Math.log(x) -0.5*x;
break;
case 13: // t-distribution (validated)
// P[0] = mu, P[1] = sigma, P[2] = df
z = 0.5*(P2 + 1);
ret = gammln(z) - gammln(0.5*P2) - Math.log(P1) -0.5*Math.log(P2*PI) -
(z)*Math.log(1 + (x-P0)*(x-P0)/(P1*P1*P2));
break;
case 14: // double exponential (laplace distribution) (validated)
ret = -Math.log(P1) - LN2 -Math.abs(x - P0)/P1;
break;
case 15: // gamma distribution (u, b)
// likelihood is just gamma likelihood
z = P0/P1;
ret = -gammln(z) -P0*Math.log(P1) + (z-1)*Math.log(x) -x/P1;
break;
case 16: // gamma (u, b, d) distribution
// EW commented out the following line to change to standard C-syntax (++) 03/07/06
//x = x - P2 + P0;
x += (-P2 + P0);
z = P0/P1;
ret = -gammln(z) -(z)*Math.log(P1) + (z-1)*Math.log(x) -x/P1;
break;
case 17: // half-normal (a, b) distribution
ret = -Math.log(P1) + 0.5*LN2 - 0.5*LNPI + -0.5*((x - P0)/P1)*((x - P0)/P1);
break;
} // end switch statement
} // end if statement
else { // discrete random numbers
switch(distID) {
case 18: // betaBinomial (a, b, N)
ret = gammln(P2+1) - gammln(x+1) - gammln(P2-x+1) + gammln(P0+x) + gammln(P1+P2-x)
- gammln(P0+P1+P2) + gammln(P0+P1) - gammln(P0) - gammln(P1);
break;
case 20:
//discrete uniform (validated)
// sample a-b, inclusive
ret = -Math.log(P1 - P0);
break;
case 21:
// poisson (validated)
ret = -P0 + x*Math.log((int)P0) -logStirling((int)x);
break;
case 22:
// negative binomial (validated)
ret = logStirling((int)(P0 + x - 1)) - logStirling((int)x) -
logStirling((int)(P0-1)) + P0*Math.log(P1) + x*Math.log(1-P1);
break;
case 23:
// binomial (validated)
ret = logStirling((int)P0) - logStirling((int)(P0 - x)) -
logStirling((int)x) + x*Math.log(P1) + (P0-x)*Math.log(1-P1);
break;
case 24:
// geometric (validated)
ret = Math.log(P0) + (x - 1)*Math.log(1-P0);
break;
} // end switch
} // end else
return ret;
}
/**
* Method gammln returns the log of the gamma function.
*
* @param xx, the value to evaluate the function at
* @return double
*/
public double gammln(double xx) {
y=x=xx;
tmp=x+5.5;
ser=1.000000000190015;
tmp -= (x+0.5)*Math.log(tmp);
for (j=0;j<=5;j++) ser += cof[j]/++y;
return (double)(-tmp+Math.log(2.5066282746310005*ser/x));
}
/**
* This method uses Stirling's approximation to generate Beta
* random variables.
*
* @param n1, the first beta parameter
* @param n2, the second beta parameter
* @return double
*/
public double betaFunction(int n1, int n2) {
// beta can be re-written as gamma(n1)*gamma(n2)/gamma(n1+n2)
return stirling(n1-1)*stirling(n2-1)/stirling(n1+n2-1);
}
/**
* This method returns Stirling's approximation to (n!).
*
* @param n, the value to find the factorial of
* @return double
*/
public double stirling(int n) {
f = 1.0;
if(n <= 10) {
// calculate the actual factorial
for(j = 1; j <= n; j++) f*=j;
return f;
}
return (ROOT2PI*Math.pow((float)n,(float)n+0.5)*Math.exp(-(double)n));
}
/**
* This function returns the log of Stirling's approximation to (n!).
*
* @param n, the integer to find the factorial of
* @return double
*/
public double logStirling(int n) {
f = 1.0;
if(n <= 10) {
// calculate the actual factorial
// f = f * j
for(j = 1; j <= n; j++) f*=j;
return Math.log(f);
}
return (LOGRT2PI + (n+0.5)*Math.log(n) -(double)n);
}
}