Example usage for org.apache.commons.math.analysis.integration TrapezoidIntegrator integrate

Introduction

In this page you can find the example usage for org.apache.commons.math.analysis.integration TrapezoidIntegrator integrate.

Prototype

public double integrate(final UnivariateRealFunction f, final double min, final double max)
        throws MaxIterationsExceededException, FunctionEvaluationException, IllegalArgumentException 

Source Link

Usage

From source file:test.dr.evomodel.substmodel.TwoStateOccupancyMarkovRewardsTest.java

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;//from www  .j  a va  2 s  .  c o m
    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");
}

---

From source file:test.dr.math.GeneralizedIntegerGammaTest.java

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 {//  w ww  .java  2  s. co  m
            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];
        assertEquals(variance, m2, tolerance * 10); // Harder to approximate
    }
}

---