import RngPack.*; /** * Write a description of class randomLibrary here. * * @author Eric Ward * @version August 31, 2006 */ public class randomLibrary { // 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; final double LN4 = 1.38629436111; final double E = 2.718281828459; public float I; public double[] JDRAW; // variables for rGas double f, v1, v2, s; // variables for poiDev double sq,alxm,g,oldm,em,t,y; // these variables are used for the gammln function private double x, 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; // these variables are used for the gamdev function double A, B, Q, D, T, C, p, u, v, w, p1, p2; // randLib variables int count; // indexing variable double rnd, min, max, gamma1, gamma2;// temporary variables static RandomElement re; static long seed; /** * Constructor for objects of class randomLibrary */ public randomLibrary() { re = new Ranmar(seed); seed=RandomSeedable.ClockSeed(); JDRAW = new double[2]; I = 1; } /** * Method randLib returns random number distribution specified. * * @param distID, the unique integer key for this distribution * @param distType, an integer describing whether this is continuous (0) * or discrete (1) * @param P0, the first parameter for this distribution * @param P1, the second parameter for this distribution * @param P2, the third parameter for this distribution * @return void */ public double randLib(int distID, int distType, double P0, double P1, double P2) { //distID is the ID of the distribution: 1 = normal, 2 = uniform, etc. //distType: continuous (0) or discrete (1) //P0 the first parameter of the distribution //P1 the second parameter of the distribution //P2 the third parameter of the distribution //*idum, the seed used for the random number generators 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 = gasdev()*P1 + P0; break; case 2:// continuous uniform (validated) // PARS[0] = low, PARS[1] = high ret = re.raw()*(P1-P0) + P0; break; case 3:// log-uniform // log-uniform: take a uniform random number, return e(r) // PARS[0] = low, PARS[1] = high ret = Math.exp(re.raw()*(P1-P0) + P0); break; case 4:// gamma distribution (validated) // generate a gamma(a,b) random variable ret = gamdev(P0,P1); break; case 5:// 3-parameter gamma distribution (validated) // generate a gamma(a,b) random variable, de-mean it // and add the location param ret = gamdev(P0,P1) - P0*P1 + P2; break; case 6:// lognormal distribution (validated) // generate a normal(0,1) random variable ret = Math.exp(gasdev()*P1 + P0); break; case 7:// exponential (validated) ret = -P0*Math.log(re.raw()); break; case 8:// inverse gamma (validated) ret = 1.0/gamdev(P0, 1.0/P1); break; case 9:// logistic (validated) ret = P0 - P1*Math.log(1.0/re.raw() - 1.0); break; case 10:// beta (validated) min = 0.0; max = 1.0; gamma1 = gamdev(1.0, P0); gamma2 = gamdev(1.0, P1); if(P0 < P1) ret = max-(max-min)*gamma2/(gamma1+gamma2); else ret = min + (max-min)*gamma1/(gamma1+gamma2); break; case 11:// f distribution (validated) //gamma1 = gamdev(2.0, P0/2.0); //gamma2 = gamdev(2.0, P1/2.0); ret = (gamdev(2.0, P0/2.0)/P0)/(gamdev(2.0, P1/2.0)/P1); break; case 12:// chi-square ret = gamdev(0.5*P0, 2.0); break; case 13: // t-distribution (validated) // P[0] = mu, P[1] = sigma, P[2] = df // draw standard normal //rnd = gasdev(); // generate a chi-square(df) RV //gamma1 = gamdev(P2/2.0, 2.0); ret = P0 + P1*gasdev()/Math.sqrt(gamdev(P2/2.0, 2.0)/P2); break; case 14: // double exponential (laplace distribution) (validated) // generate two iid (0,1) random numbers min = re.raw(); max = re.raw(); if(min > 0.5) { ret = P0 + P1*Math.log(max); } else { ret = P0 - P1*Math.log(max); } break; case 15: // gamma (u, b) distribution ret = gamdev(P0/P1,P1); break; case 16: // gamma (u, b, d) distribution // centralized gamma distribution (validated) // generate a gamma(a,b) random variable, de-mean it // and add location param ret = gamdev(0/P1,P1) - P0 + P2; break; case 17: // half-normal (a, b) distribution ret = gasdev()*P1 + P0; break; } // end switch statement } // end if statement else { // discrete random numbers switch(distID) { case 18: // beta-binomial(a,b,N) min = 0.0; max = 1.0; gamma1 = gamdev(1.0, P0); gamma2 = gamdev(1.0, P1); if(P0 < P1) rnd = max-(max-min)*gamma2/(gamma1+gamma2); else rnd = min + (max-min)*gamma1/(gamma1+gamma2); count = 0; for(j = 0; j < (int)P2; j++) { ret = re.raw(); if(ret < rnd) { count++; } } ret = count; break; case 20: //discrete uniform (validated) // sample a-b, inclusive ret = P0 + Math.ceil(re.raw()*(P1-P0+1))-1.0; break; case 21: // poisson (validated) ret = poidev(P0); break; case 22: // negative binomial (validated) ret = 0; for(j = 0; j < (int)P0; j++) { while(true) { rnd = re.raw(); if(rnd <= P1) { break; } ret++; } } break; case 23: // binomial (validated) // cycle through N times count = 0; for(j = 0; j < (int)P0; j++) { rnd = re.raw(); if(rnd < P1) { count++; } } ret = count; break; case 24: // geometric (validated) ret = 0; while(true) { ret ++; rnd = re.raw(); if(rnd < P0) { break; } } break; } // end switch } // end else return ret; } /** * Method gammln returns the log of the gamma function. * * @param xx, the value to evaluate this function at * @return value of the function */ 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)); } /** * Method getJump is used to efficiently draw normal random numbers * for the MCMC jumping process. Normally, a routine generates two * independent numbers each time it is called (and one is thrown out). * This function aims to use both draws. * * @param void * @return double JDRAW */ public double getJump() { if(I == 0) { I = 1; return JDRAW[1]; } else { I = 0; s = 2.0; do { v1 = 2*re.raw() - 1; v2 = 2*re.raw() - 1; s = v1*v1 + v2*v2; } while(s > 1); f = Math.sqrt(-2.0*Math.log(s)/s); JDRAW[0] = f*v1; JDRAW[1] = f*v2; return JDRAW[0]; } } /** * This method draws a uniform random number. * * @param void * @return double rnd */ public double unif() { return re.raw(); } /** * Method rGas draws two iid standard normal random variables, and * is taken from Numerical Recipes. * * @param z, a 2x1 array used to store the draws temporarily * @return void */ public void rGas(double z[]) { s = 2.0; do { v1 = 2*re.raw() - 1; v2 = 2*re.raw() - 1; s = v1*v1 + v2*v2; } while(s > 1); f = Math.sqrt(-2.0*Math.log(s)/s); z[0] = f*v1; z[1] = f*v2; } /** * Method gasdev draws one iid standard normal random variable. * * @param void * @return double rnd */ public double gasdev() { s = 2.0; do { v1 = 2*re.raw() - 1; v2 = 2*re.raw() - 1; s = v1*v1 + v2*v2; } while(s > 1); f = Math.sqrt(-2.0*Math.log(s)/s); return f*v1; } /** * Method poidev draws Poisson random numbers. * * @param xm, the mean value * @return double rnd */ public double poidev(double xm) { oldm = -1.0; if (xm < 12.0) { if (xm != oldm) { oldm=xm; g=Math.exp(-xm); } em = -1; t=1.0; do { ++em; t *= re.raw(); } while (t > g); } else { if (xm != oldm) { oldm=xm; sq=Math.sqrt(2.0*xm); alxm=Math.log(xm); g=xm*alxm-gammln((xm+1.0)); } do { do { y=Math.tan(PI*re.raw()); em=sq*y+xm; } while (em < 0.0); em=Math.floor(em); t=(0.9*(1.0+y*y)*Math.exp(em*alxm-gammln((em+1.0))-g)); } while (re.raw() > t); } return em; } /** * Method gamdev is the Numerical Recipes routine for creating * gamma random variables. * * @param a, normally thought of as gamma b - when a = 1, this is exponential * @param b, normally thought of as gamma a, the scale parameter * @return double rnd, */ public double gamdev(double a, double b) { // three algorithms depending on the value of a // note: a and b are switched from traditional notation A = Math.sqrt(2.0*a - 1.0); B = a - LN4; Q = a + 1.0/A; T = 4.5; D = 1.0+Math.log(T); C = 1.0 + a/(E); if(a < 1.0) { while(true) { p = re.raw(); if(p > 1.0) { y = -Math.log((C-p)/a); u = re.raw(); if(u <= Math.pow(y,a-1.0)) { return (b*y); } } else { y = Math.pow(p,1.0/a); u = re.raw(); if(u <= Math.exp(-y)) { return b*y; } } } } else if(a == 1) { return -(b)*Math.log(re.raw()); } else { while(true) { p1 = re.raw(); p2 = re.raw(); v = A*Math.log(p1/(1.0-p1)); y = (a*Math.exp(v)); z = p1*p1*p2; w = B + Q*v - y; if(w + D - T*z >= 0 || w >= Math.log(z)) { return b*y; } } } } }