Wednesday 28 November 2007

Sample from Gamma distribution in Java

Back to the Java struggles.
Today I am giving away a piece of Java code which allows one to sample a random variable from Gamma distribution.
import java.util.Random;

public class Samplers {
private static Random rng = new Random(
Calendar.getInstance().getTimeInMillis() +
Thread.currentThread().getId());

public static double sampleGamma(double k, double theta) {
boolean accept = false;
if (k < 1) {
// Weibull algorithm
double c = (1 / k);
double d = ((1 - k) * Math.pow(k, (k / (1 - k))));
double u, v, z, e, x;
do {
u = rng.nextDouble();
v = rng.nextDouble();
z = -Math.log(u);
e = -Math.log(v);
x = Math.pow(z, c);
if ((z + e) >= (d + x)) {
accept = true;
}
} while (!accept);
return (x * theta);
} else {
// Cheng's algorithm
double b = (k - Math.log(4));
double c = (k + Math.sqrt(2 * k - 1));
double lam = Math.sqrt(2 * k - 1);
double cheng = (1 + Math.log(4.5));
double u, v, x, y, z, r;
do {
u = rng.nextDouble();
v = rng.nextDouble();
y = ((1 / lam) * Math.log(v / (1 - v)));
x = (k * Math.exp(y));
z = (u * v * v);
r = (b + (c * y) - x);
if ((r >= ((4.5 * z) - cheng)) ||
(r >= Math.log(z))) {
accept = true;
}
} while (!accept);
return (x * theta);
}
}
}