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

Example 1 with UnivariateRealFunction

use of org.apache.commons.math.analysis.UnivariateRealFunction 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 UnivariateRealFunction

use of org.apache.commons.math.analysis.UnivariateRealFunction 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.");
    }
}

Example 3 with UnivariateRealFunction

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

the class GeneralizedIntegerGammaTest method testPdf.

public void testPdf() {
    for (int i = 0; i < shape1s.length; ++i) {
        final GeneralizedIntegerGammaDistribution gig = new GeneralizedIntegerGammaDistribution(shape1s[i], shape2s[i], rate1s[i], rate2s[i]);
        TrapezoidIntegrator integrator = new TrapezoidIntegrator();
        double m0 = 0.0;
        double m1 = 0.0;
        double m2 = 0.0;
        try {
            m0 = integrator.integrate(new UnivariateRealFunction() {

                @Override
                public double value(double x) throws FunctionEvaluationException {
                    final double pdf = gig.pdf(x);
                    return pdf;
                }
            }, 0.0, 1000.0);
            m1 = integrator.integrate(new UnivariateRealFunction() {

                @Override
                public double value(double x) throws FunctionEvaluationException {
                    final double pdf = gig.pdf(x);
                    return x * pdf;
                }
            }, 0.0, 1000.0);
            m2 = integrator.integrate(new UnivariateRealFunction() {

                @Override
                public double value(double x) throws FunctionEvaluationException {
                    final double pdf = gig.pdf(x);
                    return x * x * pdf;
                }
            }, 0.0, 1000.0);
        } catch (MaxIterationsExceededException e) {
            e.printStackTrace();
        } catch (FunctionEvaluationException e) {
            e.printStackTrace();
        }
        // Check normalization
        assertEquals(1.0, m0, tolerance);
        // Check mean
        double mean = shape1s[i] / rate1s[i] + shape2s[i] / rate2s[i];
        assertEquals(mean, m1, tolerance);
        // Check variance
        m2 -= m1 * m1;
        double variance = shape1s[i] / rate1s[i] / rate1s[i] + shape2s[i] / rate2s[i] / rate2s[i];
        // Harder to approximate
        assertEquals(variance, m2, tolerance * 10);
    }
}