Examples with Gamma beast.math.distributions.Gamma used on opensource projects

Example 1 with Gamma

use of beast.math.distributions.Gamma in project beast2 by CompEvol.

the class GammaTest method testPdf.

public void testPdf() throws MathException {
    final int numberOfTests = 300;
    double totErr = 0;
    double ptotErr = 0;
    int np = 0;
    double qtotErr = 0;
    Randomizer.setSeed(551);
    for (int i = 0; i < numberOfTests; i++) {
        final double mean = .01 + (3 - 0.01) * Randomizer.nextDouble();
        final double var = .01 + (3 - 0.01) * Randomizer.nextDouble();
        double scale0 = var / mean;
        double shape = mean / scale0;
        final Gamma gamma = new Gamma();
        Gamma.mode mode = Gamma.mode.values()[Randomizer.nextInt(4)];
        double other = 0;
        switch(mode) {
            case ShapeScale:
                other = scale0;
                break;
            case ShapeRate:
                other = 1 / scale0;
                break;
            case ShapeMean:
                other = scale0 * shape;
                break;
            case OneParameter:
                other = 1 / shape;
                scale0 = 1 / shape;
                break;
        }
        final double scale = scale0;
        gamma.initByName("alpha", shape + "", "beta", other + "", "mode", mode);
        final double value = Randomizer.nextGamma(shape, 1 / scale);
        final double mypdf = mypdf(value, shape, scale);
        final double pdf = gamma.density(value);
        if (Double.isInfinite(mypdf) && Double.isInfinite(pdf)) {
            continue;
        }
        assertFalse(Double.isNaN(mypdf));
        assertFalse(Double.isNaN(pdf));
        totErr += mypdf != 0 ? Math.abs((pdf - mypdf) / mypdf) : pdf;
        assertFalse("nan", Double.isNaN(totErr));
        // assertEquals("" + shape + "," + scale + "," + value, mypdf,gamma.pdf(value),1e-10);
        final double cdf = gamma.cumulativeProbability(value);
        UnivariateRealFunction f = new UnivariateRealFunction() {

            public double value(double v) throws FunctionEvaluationException {
                return mypdf(v, shape, scale);
            }
        };
        final UnivariateRealIntegrator integrator = new RombergIntegrator();
        integrator.setAbsoluteAccuracy(1e-14);
        // fail if it takes too much time
        integrator.setMaximalIterationCount(16);
        double x;
        try {
            x = integrator.integrate(f, 0, value);
            ptotErr += cdf != 0.0 ? Math.abs(x - cdf) / cdf : x;
            np += 1;
            // assertTrue("" + shape + "," + scale + "," + value + " " + Math.abs(x-cdf)/x + "> 1e-6", Math.abs(1-cdf/x) < 1e-6);
            final double q = gamma.inverseCumulativeProbability(cdf);
            qtotErr += q != 0 ? Math.abs(q - value) / q : value;
        // System.out.println(shape + ","  + scale + " " + value);
        } catch (ConvergenceException e) {
        // can't integrate , skip test
        // System.out.print(" theta(" + shape + ","  + scale + ") skipped");
        }
    // assertEquals("" + shape + "," + scale + "," + value + " " + Math.abs(q-value)/value, q, value, 1e-6);
    // System.out.print("\n" + np + ": " + mode + " " + totErr/np + " " + qtotErr/np + " " + ptotErr/np);
    }
    // System.out.println( !Double.isNaN(totErr) );
    // System.out.println(totErr);
    // bad test, but I can't find a good threshold that works for all individual cases
    assertTrue("failed " + totErr / numberOfTests, totErr / numberOfTests < 1e-7);
    assertTrue("failed " + qtotErr / numberOfTests, qtotErr / numberOfTests < 1e-10);
    assertTrue("failed " + ptotErr / np, np > 0 ? (ptotErr / np < 2e-7) : true);
}

Example 2 with Gamma

use of beast.math.distributions.Gamma in project beast2 by CompEvol.

the class SpeciesTreePrior method initAndValidate.

@Override
public void initAndValidate() {
    popFunction = popFunctionInput.get();
    popSizesBottom = popSizesBottomInput.get();
    popSizesTop = popSizesTopInput.get();
    // set up sizes of population functions
    final int speciesCount = treeInput.get().getLeafNodeCount();
    final int nodeCount = treeInput.get().getNodeCount();
    switch(popFunction) {
        case constant:
            popSizesBottom.setDimension(nodeCount);
            break;
        case linear:
            if (popSizesTop == null) {
                throw new IllegalArgumentException("topPopSize must be specified");
            }
            popSizesBottom.setDimension(speciesCount);
            popSizesTop.setDimension(nodeCount);
            break;
        case linear_with_constant_root:
            if (popSizesTop == null) {
                throw new IllegalArgumentException("topPopSize must be specified");
            }
            popSizesBottom.setDimension(speciesCount);
            popSizesTop.setDimension(nodeCount - 1);
            break;
    }
    // bottom prior = Gamma(2,Psi)
    gamma2Prior = new Gamma();
    gamma2Prior.betaInput.setValue(gammaParameterInput.get(), gamma2Prior);
    // top prior = Gamma(4,Psi)
    gamma4Prior = new Gamma();
    final RealParameter parameter = new RealParameter(new Double[] { 4.0 });
    gamma4Prior.alphaInput.setValue(parameter, gamma4Prior);
    gamma4Prior.betaInput.setValue(gammaParameterInput.get(), gamma4Prior);
    if (popFunction != TreePopSizeFunction.constant && gamma4Prior == null) {
        throw new IllegalArgumentException("Top prior must be specified when population function is not constant");
    }
// make sure the m_taxonSet is a set of taxonsets
// HACK to make Beauti initialise: skip the check here
// for (Taxon taxon : m_taxonSet.get().m_taxonset.get()) {
// if (!(taxon instanceof TaxonSet)) {
// throw new IllegalArgumentException("taxonset should be sets of taxa only, not individual taxons");
// }
// }
}

Example 3 with Gamma

use of beast.math.distributions.Gamma in project beast2 by CompEvol.

the class GammaTest method testGammaCummulative.

@Test
public void testGammaCummulative() throws Exception {
    Gamma dist = new Gamma();
    dist.initByName("alpha", "0.001", "beta", "1000.0");
    double v = dist.inverseCumulativeProbability(0.5);
    assertEquals(v, 5.244206e-299, 1e-304);
    v = dist.inverseCumulativeProbability(0.05);
    assertEquals(v, 0.0);
    v = dist.inverseCumulativeProbability(0.025);
    assertEquals(v, 0.0);
    v = dist.inverseCumulativeProbability(0.95);
    assertEquals(v, 2.973588e-20, 1e-24);
    v = dist.inverseCumulativeProbability(0.975);
    assertEquals(v, 5.679252e-09, 1e-13);
}

Example 4 with Gamma

use of beast.math.distributions.Gamma in project beast2 by CompEvol.

the class MeanOfParametricDistributionTest method testMeanOfGamma.

@Test
public void testMeanOfGamma() throws Exception {
    Gamma gamma = new Gamma();
    gamma.initByName("alpha", "100", "beta", "10");
    double mean = gamma.getMean();
    assertEquals(mean, 1000, 1e-10);
    gamma = new Gamma();
    gamma.initByName("alpha", "100", "beta", "100");
    mean = gamma.getMean();
    assertEquals(mean, 10000, 1e-10);
    gamma = new Gamma();
    gamma.initByName("alpha", "100", "beta", "10", "offset", "3");
    mean = gamma.getMean();
    assertEquals(mean, 1003, 1e-10);
}