Search in sources :

Example 1 with BetaDistribution

use of org.apache.commons.math3.distribution.BetaDistribution in project gatk by broadinstitute.

the class SliceSamplerUnitTest method testSliceSamplingOfPeakedBetaDistribution.

/**
     * Test slice sampling of a peaked 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 testSliceSamplingOfPeakedBetaDistribution() {
    rng.setSeed(RANDOM_SEED);
    final double alpha = 10.;
    final double beta = 4.;
    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);
}
Also used : BetaDistribution(org.apache.commons.math3.distribution.BetaDistribution) Mean(org.apache.commons.math3.stat.descriptive.moment.Mean) Variance(org.apache.commons.math3.stat.descriptive.moment.Variance) Test(org.testng.annotations.Test)

Example 2 with BetaDistribution

use of org.apache.commons.math3.distribution.BetaDistribution in project gatk-protected 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);
}
Also used : BetaDistribution(org.apache.commons.math3.distribution.BetaDistribution) Mean(org.apache.commons.math3.stat.descriptive.moment.Mean) Variance(org.apache.commons.math3.stat.descriptive.moment.Variance) Test(org.testng.annotations.Test)

Example 3 with BetaDistribution

use of org.apache.commons.math3.distribution.BetaDistribution in project gatk-protected by broadinstitute.

the class SliceSamplerUnitTest method testSliceSamplingOfPeakedBetaDistribution.

/**
     * Test slice sampling of a peaked 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 testSliceSamplingOfPeakedBetaDistribution() {
    rng.setSeed(RANDOM_SEED);
    final double alpha = 10.;
    final double beta = 4.;
    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);
}
Also used : BetaDistribution(org.apache.commons.math3.distribution.BetaDistribution) Mean(org.apache.commons.math3.stat.descriptive.moment.Mean) Variance(org.apache.commons.math3.stat.descriptive.moment.Variance) Test(org.testng.annotations.Test)

Example 4 with BetaDistribution

use of org.apache.commons.math3.distribution.BetaDistribution in project gatk 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);
}
Also used : BetaDistribution(org.apache.commons.math3.distribution.BetaDistribution) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Example 5 with BetaDistribution

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);
}
Also used : BetaDistribution(org.apache.commons.math3.distribution.BetaDistribution) BaseTest(org.broadinstitute.hellbender.utils.test.BaseTest) Test(org.testng.annotations.Test)

Aggregations

BetaDistribution (org.apache.commons.math3.distribution.BetaDistribution)7 Test (org.testng.annotations.Test)6 Mean (org.apache.commons.math3.stat.descriptive.moment.Mean)4 Variance (org.apache.commons.math3.stat.descriptive.moment.Variance)4 BaseTest (org.broadinstitute.hellbender.utils.test.BaseTest)2 RegDataSet (edu.neu.ccs.pyramid.dataset.RegDataSet)1