List of usage examples for org.apache.commons.math3.analysis.integration UnivariateIntegrator integrate
double integrate(int maxEval, UnivariateFunction f, double min, double max) throws TooManyEvaluationsException, MaxCountExceededException, MathIllegalArgumentException, NullArgumentException;
From source file:com.github.rinde.rinsim.scenario.generator.IntensityFunctions.java
/** * Compute the area of a sine intensity by using numerical approximation. The * range for which the area is computed is defined by a lower bound of * <code>0</code> and an upper bound of <code>length</code>. * @param s The intensity function to integrate. * @param lb The lower bound of the range. * @param ub The upper bound of the range. * @return The area./* w w w. ja v a2s . c om*/ */ public static double areaByIntegration(IntensityFunction s, double lb, double ub) { final UnivariateIntegrator ri = new RombergIntegrator(16, 32); final double val = ri.integrate(10000000, asUnivariateFunction(s), lb, ub); return val; }
From source file:gdsc.smlm.results.PeakResult.java
/** * Compute the function I1 using numerical integration. See Mortensen, et al (2010) Nature Methods 7, 377-383), SI * equation 43./*from w w w. j a v a2s . c o m*/ * * <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; }
From source file:beast.math.distribution.GammaDistributionTest.java
@Test public void testPdf() { final int numberOfTests = 100; double totErr = 0; double ptotErr = 0; int np = 0;/* www. ja v a2 s . c o m*/ 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); }
From source file:gdsc.smlm.function.PoissonGammaGaussianFunction.java
/** * Compute the likelihood/*from w ww.j a v a 2 s. co m*/ * * @param o * The observed count * @param e * The expected count * @return The likelihood */ public double likelihood(final double o, final double e) { // Use the same variables as the Python code final double cij = o; final double eta = alpha * e; // convert to photons if (sigma == 0) { // No convolution with a Gaussian. Simply evaluate for a Poisson-Gamma distribution. // Any observed count above zero if (cij > 0.0) { // The observed count converted to photons final double nij = alpha * cij; if (eta * nij > 10000) { // Approximate Bessel function i1(x) when using large x: // i1(x) ~ exp(x)/sqrt(2*pi*x) // However the entire equation is logged (creating transform), // evaluated then raised to e to prevent overflow error on // large exp(x) final double transform = 0.5 * Math.log(alpha * eta / cij) - nij - eta + 2 * Math.sqrt(eta * nij) - Math.log(twoSqrtPi * Math.pow(eta * nij, 0.25)); return FastMath.exp(transform); } else { // Second part of equation 135 return Math.sqrt(alpha * eta / cij) * FastMath.exp(-nij - eta) * Bessel.I1(2 * Math.sqrt(eta * nij)); } } else if (cij == 0.0) { return FastMath.exp(-eta); } else { return 0; } } else if (useApproximation) { return mortensenApproximation(cij, eta); } else { // This code is the full evaluation of equation 7 from the supplementary information // of the paper Chao, et al (2013) Nature Methods 10, 335-338. // It is the full evaluation of a Poisson-Gamma-Gaussian convolution PMF. final double sk = sigma; // Read noise final double g = 1.0 / alpha; // Gain final double z = o; // Observed pixel value final double vk = eta; // Expected number of photons // Compute the integral to infinity of: // exp( -((z-u)/(sqrt(2)*s)).^2 - u/g ) * sqrt(vk*u/g) .* besseli(1, 2 * sqrt(vk*u/g)) ./ u; final double vk_g = vk * alpha; // vk / g final double sqrt2sigma = Math.sqrt(2) * sk; // Specify the function to integrate UnivariateFunction f = new UnivariateFunction() { public double value(double u) { return eval(sqrt2sigma, z, vk_g, g, u); } }; // Integrate to infinity is not necessary. The convolution of the function with the // Gaussian should be adequately sampled using a nxSD around the maximum. // Find a bracket containing the maximum double lower, upper; double maxU = Math.max(1, o); double rLower = maxU; double rUpper = maxU + 1; double f1 = f.value(rLower); double f2 = f.value(rUpper); // Calculate the simple integral and the range double sum = f1 + f2; boolean searchUp = f2 > f1; if (searchUp) { while (f2 > f1) { f1 = f2; rUpper += 1; f2 = f.value(rUpper); sum += f2; } maxU = rUpper - 1; } else { // Ensure that u stays above zero while (f1 > f2 && rLower > 1) { f2 = f1; rLower -= 1; f1 = f.value(rLower); sum += f1; } maxU = (rLower > 1) ? rLower + 1 : rLower; } lower = Math.max(0, maxU - 5 * sk); upper = maxU + 5 * sk; if (useSimpleIntegration && lower > 0) { // If we are not at the zero boundary then we can use a simple integration by adding the // remaining points in the range for (double u = rLower - 1; u >= lower; u -= 1) { sum += f.value(u); } for (double u = rUpper + 1; u <= upper; u += 1) { sum += f.value(u); } } else { // Use Legendre-Gauss integrator try { final double relativeAccuracy = 1e-4; final double absoluteAccuracy = 1e-8; final int minimalIterationCount = 3; final int maximalIterationCount = 32; final int integrationPoints = 16; // Use an integrator that does not use the boundary since u=0 is undefined (divide by zero) UnivariateIntegrator i = new IterativeLegendreGaussIntegrator(integrationPoints, relativeAccuracy, absoluteAccuracy, minimalIterationCount, maximalIterationCount); sum = i.integrate(2000, f, lower, upper); } catch (TooManyEvaluationsException ex) { return mortensenApproximation(cij, eta); } } // Compute the final probability //final double f1 = z / sqrt2sigma; return (FastMath.exp(-vk) / (sqrt2pi * sk)) * (FastMath.exp(-(f1 * f1)) + sum); } }
From source file:gdsc.smlm.model.AiryPSFModel.java
private static synchronized void createAiryDistribution() { if (spline != null) return;/*from w ww . ja va2 s. c om*/ 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); }
From source file:org.apache.mahout.math.jet.random.DistributionChecks.java
public static void checkDistribution(final AbstractContinousDistribution dist, double[] x, double offset, double scale, int n) { double[] xs = Arrays.copyOf(x, x.length); for (int i = 0; i < xs.length; i++) { xs[i] = xs[i] * scale + offset;//from w ww .j a v a2s . com } Arrays.sort(xs); // collect samples double[] y = new double[n]; for (int i = 0; i < n; i++) { y[i] = dist.nextDouble(); } Arrays.sort(y); // compute probabilities for bins double[] p = new double[xs.length + 1]; double lastP = 0; for (int i = 0; i < xs.length; i++) { double thisP = dist.cdf(xs[i]); p[i] = thisP - lastP; lastP = thisP; } p[p.length - 1] = 1 - lastP; // count samples in each bin int[] k = new int[xs.length + 1]; int lastJ = 0; for (int i = 0; i < k.length - 1; i++) { int j = 0; while (j < n && y[j] < xs[i]) { j++; } k[i] = j - lastJ; lastJ = j; } k[k.length - 1] = n - lastJ; // now verify probabilities by comparing to integral of pdf UnivariateIntegrator integrator = new RombergIntegrator(); for (int i = 0; i < xs.length - 1; i++) { double delta = integrator.integrate(1000000, new UnivariateFunction() { @Override public double value(double v) { return dist.pdf(v); } }, xs[i], xs[i + 1]); Assert.assertEquals(delta, p[i + 1], 1.0e-6); } // finally compute G^2 of observed versus predicted. See http://en.wikipedia.org/wiki/G-test double sum = 0; for (int i = 0; i < k.length; i++) { if (k[i] != 0) { sum += k[i] * Math.log(k[i] / p[i] / n); } } sum *= 2; // sum is chi^2 distributed with degrees of freedom equal to number of partitions - 1 int dof = k.length - 1; // fisher's approximation is that sqrt(2*x) is approximately unit normal with mean sqrt(2*dof-1) double z = Math.sqrt(2 * sum) - Math.sqrt(2 * dof - 1); Assert.assertTrue(String.format("offset=%.3f scale=%.3f Z = %.1f", offset, scale, z), Math.abs(z) < 3); }
From source file:uk.ac.diamond.scisoft.ncd.core.DegreeOfOrientation.java
public Object[] process(Serializable buffer, Serializable axis, final int[] dimensions) { double[] parentaxis = (double[]) ConvertUtils.convert(axis, double[].class); float[] parentdata = (float[]) ConvertUtils.convert(buffer, float[].class); int size = dimensions[dimensions.length - 1]; double[] myaxis = new double[size]; double[] mydata = new double[size]; double[] cos2data = new double[size]; double[] sin2data = new double[size]; double[] sincosdata = new double[size]; for (int i = 0; i < parentaxis.length; i++) { myaxis[i] = Math.toRadians(parentaxis[i]); mydata[i] = parentdata[i];//from www .jav a 2 s. co m float cos2alpha = (float) Math.cos(2.0 * myaxis[i]); float sin2alpha = (float) Math.sin(2.0 * myaxis[i]); cos2data[i] = (1.0f + cos2alpha) * parentdata[i] / 2.0; sin2data[i] = (1.0f - cos2alpha) * parentdata[i] / 2.0; sincosdata[i] = sin2alpha * parentdata[i] / 2.0; } UnivariateInterpolator interpolator = new SplineInterpolator(); UnivariateFunction function = interpolator.interpolate(myaxis, mydata); UnivariateFunction cos2Function = interpolator.interpolate(myaxis, cos2data); UnivariateFunction sin2Function = interpolator.interpolate(myaxis, sin2data); UnivariateFunction sincosFunction = interpolator.interpolate(myaxis, sincosdata); UnivariateIntegrator integrator = new IterativeLegendreGaussIntegrator(15, BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY, BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY); try { float cos2mean = (float) integrator.integrate(INTEGRATION_POINTS, cos2Function, myaxis[0], myaxis[myaxis.length - 1]); float sin2mean = (float) integrator.integrate(INTEGRATION_POINTS, sin2Function, myaxis[0], myaxis[myaxis.length - 1]); float sincosmean = (float) integrator.integrate(INTEGRATION_POINTS, sincosFunction, myaxis[0], myaxis[myaxis.length - 1]); float norm = (float) integrator.integrate(INTEGRATION_POINTS, function, myaxis[0], myaxis[myaxis.length - 1]); cos2mean /= norm; sin2mean /= norm; sincosmean /= norm; float result = (float) Math.sqrt(Math.pow(cos2mean - sin2mean, 2) - 4.0 * sincosmean * sincosmean); double angle = MathUtils.normalizeAngle(Math.atan2(2.0 * sincosmean, cos2mean - sin2mean) / 2.0, Math.PI); Object[] output = new Object[] { new float[] { result }, new float[] { (float) Math.toDegrees(angle) }, new float[] { (float) (result * Math.cos(angle)), (float) (result * Math.sin(angle)) }, }; return output; } catch (TooManyEvaluationsException e) { return new Object[] { new float[] { Float.NaN }, new double[] { Double.NaN } }; } catch (MaxCountExceededException e) { return new Object[] { new float[] { Float.NaN }, new double[] { Double.NaN } }; } }
From source file:uk.ac.diamond.scisoft.ncd.core.SaxsInvariant.java
public Object[] process(Serializable buffer, Serializable errors, Serializable axis, final int[] dimensions) { double[] parentaxis = (double[]) ConvertUtils.convert(axis, double[].class); float[] parentdata = (float[]) ConvertUtils.convert(buffer, float[].class); double[] parenterrors = (double[]) ConvertUtils.convert(errors, double[].class); int shift = (parentaxis[0] > 0 ? 1 : 0); int size = dimensions[dimensions.length - 1] + shift; double[] myaxis = new double[size]; double[] mydata = new double[size]; double[] myerrors = new double[size]; if (shift > 0) { myaxis[0] = 0.0;// w ww. j a va2 s. c o m mydata[0] = 0.0; myerrors[0] = 0.0; } for (int i = 0; i < parentaxis.length; i++) { myaxis[i + shift] = parentaxis[i]; mydata[i + shift] = parentdata[i] * parentaxis[i] * parentaxis[i]; myerrors[i + shift] = parenterrors[i] * Math.pow(parentaxis[i], 4); } UnivariateInterpolator interpolator = new SplineInterpolator(); UnivariateFunction function = interpolator.interpolate(myaxis, mydata); UnivariateIntegrator integrator = new IterativeLegendreGaussIntegrator(15, BaseAbstractUnivariateIntegrator.DEFAULT_RELATIVE_ACCURACY, BaseAbstractUnivariateIntegrator.DEFAULT_ABSOLUTE_ACCURACY); try { float result = (float) integrator.integrate(INTEGRATION_POINTS, function, 0.0, myaxis[myaxis.length - 1]); IDataset data = new FloatDataset(parentdata, dimensions); IDataset qaxis = new DoubleDataset(parentaxis, dimensions); PorodPlotData porodPlotData = (PorodPlotData) SaxsAnalysisPlotType.POROD_PLOT.getSaxsPlotDataObject(); SimpleRegression regression = porodPlotData.getPorodPlotParameters(data.squeeze(), qaxis.squeeze()); Amount<Dimensionless> c4 = porodPlotData.getC4(regression); result += (float) (c4.getEstimatedValue() / myaxis[myaxis.length - 1]); double error = 0.0; for (int i = 0; i < myaxis.length; i++) { int idx1 = Math.max(0, i - 1); int idx2 = Math.min(myaxis.length - 1, i + 1); error += Math.pow((myaxis[idx2] - myaxis[idx1]), 2) * myerrors[i] / 4.0; } error += Math.pow(c4.getAbsoluteError() / myaxis[myaxis.length - 1], 2); return new Object[] { new float[] { result }, new double[] { error } }; } catch (TooManyEvaluationsException e) { return new Object[] { new float[] { Float.NaN }, new double[] { Double.NaN } }; } catch (MaxCountExceededException e) { return new Object[] { new float[] { Float.NaN }, new double[] { Double.NaN } }; } }