/** * 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); } }