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

Introduction

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

Prototype

public TrapezoidIntegrator() 

Source Link

Document

Construct an integrator.

Usage

From source file:sphericalGeo.SphericalDiffusionModel.java

/** anglerange[0] = lower bound, anglerange[1] = upper bound of range for angle2
 * **///ww  w  .j  a  v  a  2 s .  co m
public double[] sample(double[] start, double time, double precision, double[] angleRange)
        throws ConvergenceException, FunctionEvaluationException, IllegalArgumentException {

    // first, sample an angle from the spherical diffusion density
    final double inverseVariance = precision / time;

    UnivariateRealFunction function = new UnivariateRealFunction() {
        @Override
        public double value(double x) throws FunctionEvaluationException {
            double logR = -x * x * inverseVariance / 2.0 + 0.5 * Math.log(x * Math.sin(x) * inverseVariance);

            return Math.exp(logR);
        }
    };

    UnivariateRealIntegrator integrator = new TrapezoidIntegrator();
    integrator.setAbsoluteAccuracy(1.0e-10);
    integrator.setRelativeAccuracy(1.0e-14);
    integrator.setMinimalIterationCount(2);
    integrator.setMaximalIterationCount(25);

    double[] cumulative = new double[NR_OF_STEPS];
    for (int i = 1; i < NR_OF_STEPS; i++) {
        cumulative[i] = cumulative[i - 1]
                + integrator.integrate(function, (i - 1) * Math.PI / NR_OF_STEPS, i * Math.PI / NR_OF_STEPS);
    }

    // normalise 
    double sum = cumulative[NR_OF_STEPS - 1];
    for (int i = 0; i < NR_OF_STEPS; i++) {
        cumulative[i] /= sum;
    }

    int i = Randomizer.randomChoice(cumulative);
    double angle = i * Math.PI / NR_OF_STEPS;

    // now we have an angle, use this to rotate the point [0,0] over
    // this angle in a random direction angle2
    double angle2 = angleRange[0] + Randomizer.nextDouble() * (angleRange[1] - angleRange[0]);
    //angleRange[0] = angle2;

    double[] xC = new double[] { Math.cos(angle), Math.sin(angle) * Math.cos(angle2),
            Math.sin(angle) * Math.sin(angle2) };

    // convert back to latitude, longitude relative to (lat=0, long=0)
    double[] xL = cartesian2Sperical(xC, true);
    //double [] sC = spherical2Cartesian(start[0], start[1]);
    double[] position = reverseMap(xL[0], xL[1], start[0], start[1]);
    return position;
}

---

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;//w ww . j  av  a2  s  . com
    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 {/*from w w w . ja v  a  2  s .c  o 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
    }
}

---