package edu.stanford.nlp.stats;
import java.util.Random;
/**
* Represents a Gamma distribution. The way that samples are drawn is
* stolen from Yee Whye Teh's code. It won't give the probability of a variable because
* gamma is a continuous distribution. should it give the mass at that point?
*
* @author Jenny Finkel
*/
public class Gamma implements ProbabilityDistribution<Double> {
/**
*
*/
private static final long serialVersionUID = -2992079318379176178L;
public final double alpha;
public Gamma(double alpha) {
this.alpha = alpha;
}
public Double drawSample(Random random) {
return drawSample(random, alpha);
}
public static Double drawSample(Random random, double alpha) {
if (alpha <= 0.0) {
/* Not well defined, set to zero and skip. */
return 0.0;
} else if ( alpha == 1.0 ) {
/* Exponential */
return -Math.log(Math.random());
} else if (alpha < 1.0) {
/* Use Johnks generator */
double cc = 1.0 / alpha;
double dd = 1.0 / (1.0-alpha);
while (true) {
double xx = Math.pow(Math.random(), cc);
double yy = xx + Math.pow(Math.random(), dd);
if (yy <= 1.0) {
return -Math.log(Math.random()) * xx / yy;
}
}
} else {
/* Use bests algorithm */
double bb = alpha - 1.0;
double cc = 3.0 * alpha - 0.75;
while (true) {
double uu = Math.random();
double vv = Math.random();
double ww = uu * (1.0 - uu);
double yy = Math.sqrt(cc / ww) * (uu - 0.5);
double xx = bb + yy;
if (xx >= 0) {
double zz = 64.0 * ww * ww * ww * vv * vv;
if ( ( zz <= (1.0 - 2.0 * yy * yy / xx) ) ||
( Math.log(zz) <= 2.0 * (bb * Math.log(xx / bb) - yy) ) ) {
return xx;
}
}
}
}
}
// Generalized Random sampling.
public static double drawSample(Random r, double a, double b) {
return drawSample(r,a)*b;
}
public double probabilityOf(Double x) {
return 0.0; // cos its not discrete
}
public double logProbabilityOf(Double x) {
return 0.0; // cos its not discrete
}
@Override
public int hashCode() {
return (new Double(alpha)).hashCode();
}
@Override
public boolean equals(Object o) {
if (!(o instanceof Gamma)) { return false; }
return ((Gamma)o).alpha == alpha;
}
}