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;
}
}
}
}
}