import org.apache.commons.math3.analysis.integration.UnivariateIntegrator; //导入依赖的package包/类
private void cumulativeProbabilityIsOneWithRealData(final double mu, double min, double max, boolean test)
{
	double p = 0;

	UnivariateIntegrator in = new SimpsonIntegrator();

	p = in.integrate(20000, new UnivariateFunction()
	{
		public double value(double x)
		{
			double v;
			v = PoissonCalculator.likelihood(mu, x);
			//v = pgf.probability(x, mu);
			//System.out.printf("x=%f, v=%f\n", x, v);
			return v;
		}
	}, min, max);

	System.out.printf("mu=%f, p=%f\n", mu, p);
	if (test)
	{
		Assert.assertEquals(String.format("mu=%f", mu), P_LIMIT, p, 0.02);
	}
}
 

import org.apache.commons.math3.analysis.integration.UnivariateIntegrator; //导入依赖的package包/类
public double getVelocity(double startSpeed, double startTime, double endTime){
    UnivariateIntegrator integrator = new SimpsonIntegrator();
    double integralAddition = 0;
    if (!dEqual(startTime,endTime)) {
        try {
            integralAddition = startTime * integrator.integrate(10, v -> getAccel(v), startTime, endTime);
        } catch (Exception e) {
            System.err.println(e.getMessage());
        }
    }
    return startSpeed + integralAddition;
}
 

import org.apache.commons.math3.analysis.integration.UnivariateIntegrator; //导入依赖的package包/类
private void cumulativeProbabilityIsOneWithRealData(final double mu, double min, double max, boolean test)
{
	double p = 0;

	// Test using a standard Poisson-Gaussian convolution
	//min = -max;
	//final PoissonGaussianFunction pgf = PoissonGaussianFunction.createWithVariance(1, 1, VAR);

	UnivariateIntegrator in = new SimpsonIntegrator();

	p = in.integrate(20000, new UnivariateFunction()
	{
		public double value(double x)
		{
			double v;
			v = SCMOSLikelihoodWrapper.likelihood(mu, VAR, G, O, x);
			//v = pgf.probability(x, mu);
			//System.out.printf("x=%f, v=%f\n", x, v);
			return v;
		}
	}, min, max);

	//System.out.printf("mu=%f, p=%f\n", mu, p);
	if (test)
	{
		Assert.assertEquals(String.format("mu=%f", mu), P_LIMIT, p, 0.02);
	}
}
 

import org.apache.commons.math3.analysis.integration.UnivariateIntegrator; //导入依赖的package包/类
private static synchronized void createAiryDistribution()
{
	if (spline != null)
		return;
	
	final double relativeAccuracy = 1e-4;
	final double absoluteAccuracy = 1e-8;
	final int minimalIterationCount = 3;
	final int maximalIterationCount = 32;

	UnivariateIntegrator integrator = new SimpsonIntegrator(relativeAccuracy, absoluteAccuracy,
			minimalIterationCount, maximalIterationCount);
	UnivariateFunction f = new UnivariateFunction()
	{
		public double value(double x)
		{
			// The pattern profile is in one dimension. 
			// Multiply by the perimeter of a circle to convert to 2D volume then normalise by 4 pi
			//return AiryPattern.intensity(x) * 2 * Math.PI * x / (4 * Math.PI);
			return AiryPattern.intensity(x) * 0.5 * x;
		}
	};

	// Integrate up to a set number of dark rings
	int samples = 1000;
	final double step = RINGS[SAMPLE_RINGS] / samples;
	double to = 0, from = 0;
	r = new double[samples + 1];
	sum = new double[samples + 1];
	for (int i = 1; i < sum.length; i++)
	{
		from = to;
		r[i] = to = step * i;
		sum[i] = integrator.integrate(2000, f, from, to) + sum[i - 1];
	}

	if (DoubleEquality.relativeError(sum[samples], POWER[SAMPLE_RINGS]) > 1e-3)
		throw new RuntimeException("Failed to create the Airy distribution");

	SplineInterpolator si = new SplineInterpolator();
	spline = si.interpolate(sum, r);
}
 

import org.apache.commons.math3.analysis.integration.UnivariateIntegrator; //导入依赖的package包/类
/**
 * Compute the function I1 using numerical integration. See Mortensen, et al (2010) Nature Methods 7, 377-383), SI
 * equation 43.
 * 
 * <pre>
 * I1 = 1 + sum [ ln(t) / (1 + t/rho) ] dt
 *    = - sum [ t * ln(t) / (t + rho) ] dt
 * </pre>
 * 
 * Where sum is the integral between 0 and 1. In the case of rho=0 the function returns 1;
 * 
 * @param rho
 * @param integrationPoints
 *            the number of integration points for the LegendreGaussIntegrator
 * @return the I1 value
 */
private static double computeI1(final double rho, int integrationPoints)
{
	if (rho == 0)
		return 1;

	final double relativeAccuracy = 1e-4;
	final double absoluteAccuracy = 1e-8;
	final int minimalIterationCount = 3;
	final int maximalIterationCount = 32;

	// Use an integrator that does not use the boundary since log(0) is undefined.
	UnivariateIntegrator i = new IterativeLegendreGaussIntegrator(integrationPoints, relativeAccuracy,
			absoluteAccuracy, minimalIterationCount, maximalIterationCount);

	// Specify the function to integrate
	UnivariateFunction f = new UnivariateFunction()
	{
		public double value(double x)
		{
			return x * Math.log(x) / (x + rho);
		}
	};
	final double i1 = -i.integrate(2000, f, 0, 1);
	//System.out.printf("I1 = %f (%d)\n", i1, i.getEvaluations());

	// The function requires more evaluations and sometimes does not converge,
	// presumably because log(x) significantly changes as x -> 0 where as x log(x) in the function above 
	// is more stable

	//		UnivariateFunction f2 = new UnivariateFunction()
	//		{
	//			@Override
	//			public double value(double x)
	//			{
	//				return Math.log(x) / ( 1 + x / rho);
	//			}
	//		};
	//		double i2 = 1 + i.integrate(2000, f2, 0, 1);
	//		System.out.printf("I1 (B) = %f (%d)\n", i2, i.getEvaluations());

	return i1;
}
 

import org.apache.commons.math3.analysis.integration.UnivariateIntegrator; //导入依赖的package包/类
@Test
public void testPdf() {

       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);
           UnivariateFunction f = new UnivariateFunction() {
               public double value(double v) {
                   return mypdf(v, shape, scale);
               }
           };
           final UnivariateIntegrator integrator = new RombergIntegrator(MachineAccuracy.SQRT_EPSILON, 1e-14, 1, 16);

           double x;
           try {
               x = integrator.integrate(16, f, 0.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(MaxCountExceededException 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);
}