List of usage examples for org.apache.commons.math3.analysis.integration IterativeLegendreGaussIntegrator IterativeLegendreGaussIntegrator
public IterativeLegendreGaussIntegrator(final int n, final double relativeAccuracy, final double absoluteAccuracy, final int minimalIterationCount, final int maximalIterationCount) throws NotStrictlyPositiveException, NumberIsTooSmallException
From source file:gdsc.smlm.function.PoissonGammaGaussianFunction.java
/** * Compute the likelihood//from www . ja va2s . com * * @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.results.PeakResult.java
/** * Compute the function I1 using numerical integration. See Mortensen, et al (2010) Nature Methods 7, 377-383), SI * equation 43.//w w w .j av a 2s.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; }