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

Example 1 with TrapezoidIntegrator

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

the class TwoStateOccupancyMarkovRewardsTest method run.

//    public void testTwoStateSericolaRewards1() {
//        final double rate = 0.0015;
////        final double prop = 0.5;
//        final double prop = 0.66666;
//
//        final double branchLength = 2000.0;
//        final boolean print = false;
//
////        MarkovReward markovReward = createMarkovReward(rate, prop);
//        MarkovReward markovReward = createSericolaMarkovReward(rate, prop);
//
//        run(markovReward, rate, prop, branchLength, print, 1000);
//    }
//    public void testTwoStateSericolaRewards2() {
//        final double rate = 0.0015;
//        final double prop = 0.5;
////        final double prop = 0.66666;
//        final double branchLength = 1000.0;
//        final boolean print = false;
//
//        MarkovReward markovReward = createMarkovReward(rate, prop);
////        MarkovReward markovReward = createSericolaMarkovReward(rate, prop);
//
//        run(markovReward, rate, prop, branchLength, print, 1000);
//    }
//    public void testLatentStateBranchRateModel() throws FunctionEvaluationException, MaxIterationsExceededException {
//
//        LatentStateBranchRateModel model = new LatentStateBranchRateModel(
//                new Parameter.Default(0.001), new Parameter.Default(0.5));
//
//        TrapezoidIntegrator integator = new TrapezoidIntegrator();
//
//        final double branchLength = 2000;
//        double integral = integator.integrate(new LatentStateDensityFunction(model, branchLength), 0.0, 1.0);
//
//        System.out.println("testLatentStateBeanchRateModel");
//        System.out.println("Integral = " + integral);
//
//        assertEquals(integral, 1.0, tolerance);
//    }
private void run(MarkovReward markovReward, double rate, double prop, double branchLength, boolean print, int length) {
    DensityFunction densityFunction = new DensityFunction(markovReward, branchLength, rate, prop);
    final double step = branchLength / length;
    int i = 0;
    double sum = 0.0;
    double modeY = 0.0;
    double modeX = 0.0;
    for (double x = 0.0; x <= branchLength; x += step, ++i) {
        double density = 0;
        density = densityFunction.value(x);
        if (x == 0.0) {
            modeY = density;
        } else {
            if (density > modeY) {
                modeY = density;
                modeX = x;
            }
        }
        if (x == 0.0 || x == branchLength) {
            sum += density;
        } else {
            sum += 2.0 * density;
        }
        if (print) {
            System.out.println(i + "\t" + String.format("%3.2f", x) + "\t" + String.format("%5.3e", density));
        }
    }
    sum *= (branchLength / 2.0 / length);
    // TODO Normalization is missing in LatentBranchRateModel
    System.out.println("branchLength = " + branchLength);
    System.out.println("rate = " + rate);
    System.out.println("prop = " + prop);
    System.out.println("Integral = " + sum);
    System.out.println("Mode = " + String.format("%3.2e", modeY) + " at " + modeX);
    assertEquals(sum, 1.0, tolerance);
    TrapezoidIntegrator integrator = new TrapezoidIntegrator();
    double integral = 0.0;
    try {
        integral = integrator.integrate(new UnitDensityFunction(markovReward, branchLength, rate, prop), 0.0, 1.0);
    } catch (MaxIterationsExceededException e) {
        e.printStackTrace();
    } catch (FunctionEvaluationException e) {
        e.printStackTrace();
    }
    System.out.println("unt int = " + integral);
    assertEquals(integral, 1.0, tolerance);
    System.out.println("\n");
}

Example 2 with TrapezoidIntegrator

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