[beast-mcmc] r6739 committed - Test code for gamma distribution.

[beast...@googlecode.com]

Revision: 6739
Author: bael...@gmail.com
Date: Thu Apr 16 14:36:12 2015 UTC
Log: Test code for gamma distribution.
https://code.google.com/p/beast-mcmc/source/detail?r=6739

Modified:
/trunk/src/dr/math/distributions/GammaDistribution.java

=======================================
--- /trunk/src/dr/math/distributions/GammaDistribution.java Fri Feb 15
19:58:20 2013 UTC
+++ /trunk/src/dr/math/distributions/GammaDistribution.java Thu Apr 16
14:36:12 2015 UTC
@@ -62,7 +62,8 @@

if (TRY_COLT) {
randomEngine = new MersenneTwister(MathUtils.nextInt());
- coltGamma = new Gamma(shape, scale, randomEngine);
+ System.out.println("Colt Gamma(" + shape + "," + scale + ")");
+ coltGamma = new Gamma(shape, 1.0/scale, randomEngine);
}
}

@@ -83,15 +84,27 @@
}

public double pdf(double x) {
- return pdf(x, shape, scale);
+ if (TRY_COLT) {
+ return coltGamma.pdf(x);
+ } else {
+ return pdf(x, shape, scale);
+ }
}

public double logPdf(double x) {
- return logPdf(x, shape, scale);
+ if (TRY_COLT) {
+ return Math.log(coltGamma.pdf(x));
+ } else {
+ return logPdf(x, shape, scale);
+ }
}

public double cdf(double x) {
- return cdf(x, shape, scale);
+ if (TRY_COLT) {
+ return coltGamma.cdf(x);
+ } else {
+ return cdf(x, shape, scale);
+ }
}

public double quantile(double y) {
@@ -535,6 +548,22 @@
System.out.println(e.getMessage());
}
}
+
+ GammaDistribution gamma = new GammaDistribution(0.01,100.0);
+ double[] samples = new double[100000];
+ double sum = 0.0;
+ for (int i = 0; i < samples.length; i++) {
+ samples[i] = gamma.nextGamma();
+ sum += samples[i];
+ }
+ double mean = sum/(double)samples.length;
+ System.out.println("Mean = " + mean);
+ double variance = 0.0;
+ for (int i = 0; i < samples.length; i++) {
+ variance += Math.pow((samples[i] - mean),2.0);
+ }
+ variance = variance/(double)samples.length;
+ System.out.println("Variance = " + variance);

// System.out
// .println("K-S critical values: 1.22(10%), 1.36(5%),
1.63(1%)\n");