use of org.apache.commons.math3.distribution.BetaDistribution in project gatk by broadinstitute.
the class SliceSamplerUnitTest method testSliceSamplingOfMonotonicBetaDistribution.
/**
* Test slice sampling of a monotonic beta distribution as an example of sampling of a bounded random variable.
* Checks that input mean and variance are recovered by 10000 samples to a relative error of 0.5% and 2%,
* respectively.
*/
@Test
public void testSliceSamplingOfMonotonicBetaDistribution() {
rng.setSeed(RANDOM_SEED);
final double alpha = 10.;
final double beta = 1.;
final BetaDistribution betaDistribution = new BetaDistribution(alpha, beta);
final Function<Double, Double> betaLogPDF = betaDistribution::logDensity;
final double xInitial = 0.5;
final double xMin = 0.;
final double xMax = 1.;
final double width = 0.1;
final int numSamples = 10000;
final SliceSampler betaSampler = new SliceSampler(rng, betaLogPDF, xMin, xMax, width);
final double[] samples = Doubles.toArray(betaSampler.sample(xInitial, numSamples));
final double mean = betaDistribution.getNumericalMean();
final double variance = betaDistribution.getNumericalVariance();
final double sampleMean = new Mean().evaluate(samples);
final double sampleVariance = new Variance().evaluate(samples);
Assert.assertEquals(relativeError(sampleMean, mean), 0., 0.005);
Assert.assertEquals(relativeError(sampleVariance, variance), 0., 0.02);
}
use of org.apache.commons.math3.distribution.BetaDistribution in project gatk-protected by broadinstitute.
the class SomaticLikelihoodsEngineUnitTest method testDirichletNormalization.
@Test
public void testDirichletNormalization() {
// a Dirichlet with two parameters is a Beta
final double a = 11;
final double b = 36;
final double[] params1 = new double[] { a, b };
final double normalization = Math.pow(10, SomaticLikelihoodsEngine.log10DirichletNormalization(params1));
final BetaDistribution bd = new BetaDistribution(a, b);
for (final double x : new double[] { 0.1, 0.3, 0.6 }) {
Assert.assertEquals(bd.density(x), normalization * Math.pow(x, a - 1) * Math.pow(1 - x, b - 1), 1e-6);
}
// a Dirichlet with parameters equal to 1 is flat, so the normalization is 1/the hypervolume of the simplex
// which is d! where d is the dimension
final double[] params2 = new double[] { 1, 1, 1, 1 };
final double normalization2 = Math.pow(10, SomaticLikelihoodsEngine.log10DirichletNormalization(params2));
Assert.assertEquals(normalization2, 6, 1e-6);
}
Aggregations