Examples with UnivariateRealIntegrator org.apache.commons.math.analysis.integration.UnivariateRealIntegrator used on opensource projects

Example 1 with UnivariateRealIntegrator

use of org.apache.commons.math.analysis.integration.UnivariateRealIntegrator in project beast-mcmc by beast-dev.

the class GammaDistributionTest method testPdf.

public void testPdf() throws FunctionEvaluationException {
    final int numberOfTests = 100;
    double totErr = 0;
    double ptotErr = 0;
    int np = 0;
    double qtotErr = 0;
    Random random = new Random(37);
    for (int i = 0; i < numberOfTests; i++) {
        final double mean = .01 + (3 - 0.01) * random.nextDouble();
        final double var = .01 + (3 - 0.01) * random.nextDouble();
        final double scale = var / mean;
        final double shape = mean / scale;
        final GammaDistribution gamma = new GammaDistribution(shape, scale);
        final double value = gamma.nextGamma();
        final double mypdf = mypdf(value, shape, scale);
        final double pdf = gamma.pdf(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.cdf(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);
        //System.out.println(shape + ","  + scale + " " + value);
        } catch (ConvergenceException e) {
        // can't integrate , skip test
        //  System.out.println(shape + ","  + scale + " skipped");
        }
        final double q = gamma.quantile(cdf);
        qtotErr += q != 0 ? Math.abs(q - value) / q : value;
    // assertEquals("" + shape + "," + scale + "," + value + " " + Math.abs(q-value)/value, q, value, 1e-6);
    }
    //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 " + ptotErr / np, np > 0 ? (ptotErr / np < 1e-5) : true);
    assertTrue("failed " + qtotErr / numberOfTests, qtotErr / numberOfTests < 1e-7);
}

Example 2 with UnivariateRealIntegrator

use of org.apache.commons.math.analysis.integration.UnivariateRealIntegrator in project beast-mcmc by beast-dev.

the class ReflectedNormalDistribution method main.

public static void main(String[] args) {
    //        final ReflectedNormalDistribution rnd = new ReflectedNormalDistribution(2, 2, 0.5, 2, 1e-6);
    final ReflectedNormalDistribution rnd = new ReflectedNormalDistribution(2, 2, 1, 2, 1e-6);
    rnd.pdf(1);
    UnivariateRealFunction f = new UnivariateRealFunction() {

        public double value(double v) throws FunctionEvaluationException {
            return rnd.pdf(v);
        }
    };
    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, rnd.lower, rnd.upper);
        //                      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);
        System.out.println("Integrated pdf = " + x);
    //System.out.println(shape + ","  + scale + " " + value);
    } catch (ConvergenceException e) {
    // can't integrate , skip test
    //  System.out.println(shape + ","  + scale + " skipped");
    } catch (FunctionEvaluationException e) {
        throw new RuntimeException("I have no idea why I am here.");
    }
}