use of org.apache.commons.math3.analysis.integration.IterativeLegendreGaussIntegrator in project GDSC-SMLM by aherbert.
the class PoissonGammaGaussianFunction method likelihood.
/**
* Compute the likelihood
*
* @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;
// convert to photons
final double eta = alpha * e;
if (sigma == 0) {
// No convolution with a Gaussian. Simply evaluate for a Poisson-Gamma distribution.
final double p;
// Any observed count above zero
if (cij > 0.0) {
// The observed count converted to photons
final double nij = alpha * cij;
// The limit on eta * nij is therefore (709/2)^2 = 125670.25
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));
p = FastMath.exp(transform);
} else {
// Second part of equation 135
p = Math.sqrt(alpha * eta / cij) * FastMath.exp(-nij - eta) * Bessel.I1(2 * Math.sqrt(eta * nij));
}
} else if (cij == 0.0) {
p = FastMath.exp(-eta);
} else {
p = 0;
}
return (p > minimumProbability) ? p : minimumProbability;
} 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.
// Read noise
final double sk = sigma;
// Gain
final double g = 1.0 / alpha;
// Observed pixel value
final double z = o;
// Expected number of photons
final double vk = eta;
// 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;
// vk / g
final double vk_g = vk * alpha;
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) {
// 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;
final double p = (FastMath.exp(-vk) / (sqrt2pi * sk)) * (FastMath.exp(-(f1 * f1)) + sum);
return (p > minimumProbability) ? p : minimumProbability;
}
}
use of org.apache.commons.math3.analysis.integration.IterativeLegendreGaussIntegrator in project GDSC-SMLM by aherbert.
the class PeakResult method computeI1.
/**
* 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);
return i1;
}
Aggregations