Examples with Gaussian org.apache.commons.math3.analysis.function.Gaussian used on opensource projects

Example 1 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class EMGainAnalysis method fit.

/**
	 * Fit the EM-gain distribution (Gaussian * Gamma)
	 * 
	 * @param h
	 *            The distribution
	 */
private void fit(int[] h) {
    final int[] limits = limits(h);
    final double[] x = getX(limits);
    final double[] y = getY(h, limits);
    Plot2 plot = new Plot2(TITLE, "ADU", "Frequency");
    double yMax = Maths.max(y);
    plot.setLimits(limits[0], limits[1], 0, yMax);
    plot.setColor(Color.black);
    plot.addPoints(x, y, Plot2.DOT);
    Utils.display(TITLE, plot);
    // Estimate remaining parameters. 
    // Assuming a gamma_distribution(shape,scale) then mean = shape * scale
    // scale = gain
    // shape = Photons = mean / gain
    double mean = getMean(h) - bias;
    // Note: if the bias is too high then the mean will be negative. Just move the bias.
    while (mean < 0) {
        bias -= 1;
        mean += 1;
    }
    double photons = mean / gain;
    if (simulate)
        Utils.log("Simulated bias=%d, gain=%s, noise=%s, photons=%s", (int) _bias, Utils.rounded(_gain), Utils.rounded(_noise), Utils.rounded(_photons));
    Utils.log("Estimate bias=%d, gain=%s, noise=%s, photons=%s", (int) bias, Utils.rounded(gain), Utils.rounded(noise), Utils.rounded(photons));
    final int max = (int) x[x.length - 1];
    double[] g = pdf(max, photons, gain, noise, (int) bias);
    plot.setColor(Color.blue);
    plot.addPoints(x, g, Plot2.LINE);
    Utils.display(TITLE, plot);
    // Perform a fit
    CustomPowellOptimizer o = new CustomPowellOptimizer(1e-6, 1e-16, 1e-6, 1e-16);
    double[] startPoint = new double[] { photons, gain, noise, bias };
    int maxEval = 3000;
    String[] paramNames = { "Photons", "Gain", "Noise", "Bias" };
    // Set bounds
    double[] lower = new double[] { 0, 0.5 * gain, 0, bias - noise };
    double[] upper = new double[] { 2 * photons, 2 * gain, gain, bias + noise };
    // Restart until converged.
    // TODO - Maybe fix this with a better optimiser. This needs to be tested on real data.
    PointValuePair solution = null;
    for (int iter = 0; iter < 3; iter++) {
        IJ.showStatus("Fitting histogram ... Iteration " + iter);
        try {
            // Basic Powell optimiser
            MultivariateFunction fun = getFunction(limits, y, max, maxEval);
            PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(fun), GoalType.MINIMIZE, new InitialGuess((solution == null) ? startPoint : solution.getPointRef()));
            if (solution == null || optimum.getValue() < solution.getValue()) {
                double[] point = optimum.getPointRef();
                // Check the bounds
                for (int i = 0; i < point.length; i++) {
                    if (point[i] < lower[i] || point[i] > upper[i]) {
                        throw new RuntimeException(String.format("Fit out of of estimated range: %s %f", paramNames[i], point[i]));
                    }
                }
                solution = optimum;
            }
        } catch (Exception e) {
            IJ.log("Powell error: " + e.getMessage());
            if (e instanceof TooManyEvaluationsException) {
                maxEval = (int) (maxEval * 1.5);
            }
        }
        try {
            // Bounded Powell optimiser
            MultivariateFunction fun = getFunction(limits, y, max, maxEval);
            MultivariateFunctionMappingAdapter adapter = new MultivariateFunctionMappingAdapter(fun, lower, upper);
            PointValuePair optimum = o.optimize(new MaxEval(maxEval), new ObjectiveFunction(adapter), GoalType.MINIMIZE, new InitialGuess(adapter.boundedToUnbounded((solution == null) ? startPoint : solution.getPointRef())));
            double[] point = adapter.unboundedToBounded(optimum.getPointRef());
            optimum = new PointValuePair(point, optimum.getValue());
            if (solution == null || optimum.getValue() < solution.getValue()) {
                solution = optimum;
            }
        } catch (Exception e) {
            IJ.log("Bounded Powell error: " + e.getMessage());
            if (e instanceof TooManyEvaluationsException) {
                maxEval = (int) (maxEval * 1.5);
            }
        }
    }
    IJ.showStatus("");
    IJ.showProgress(1);
    if (solution == null) {
        Utils.log("Failed to fit the distribution");
        return;
    }
    double[] point = solution.getPointRef();
    photons = point[0];
    gain = point[1];
    noise = point[2];
    bias = (int) Math.round(point[3]);
    String label = String.format("Fitted bias=%d, gain=%s, noise=%s, photons=%s", (int) bias, Utils.rounded(gain), Utils.rounded(noise), Utils.rounded(photons));
    Utils.log(label);
    if (simulate) {
        Utils.log("Relative Error bias=%s, gain=%s, noise=%s, photons=%s", Utils.rounded(relativeError(bias, _bias)), Utils.rounded(relativeError(gain, _gain)), Utils.rounded(relativeError(noise, _noise)), Utils.rounded(relativeError(photons, _photons)));
    }
    // Show the PoissonGammaGaussian approximation
    double[] f = null;
    if (showApproximation) {
        f = new double[x.length];
        PoissonGammaGaussianFunction fun = new PoissonGammaGaussianFunction(1.0 / gain, noise);
        final double expected = photons * gain;
        for (int i = 0; i < f.length; i++) {
            f[i] = fun.likelihood(x[i] - bias, expected);
        //System.out.printf("x=%d, g=%f, f=%f, error=%f\n", (int) x[i], g[i], f[i],
        //		gdsc.smlm.fitting.utils.DoubleEquality.relativeError(g[i], f[i]));
        }
        yMax = Maths.maxDefault(yMax, f);
    }
    // Replot
    g = pdf(max, photons, gain, noise, (int) bias);
    plot = new Plot2(TITLE, "ADU", "Frequency");
    plot.setLimits(limits[0], limits[1], 0, yMax * 1.05);
    plot.setColor(Color.black);
    plot.addPoints(x, y, Plot2.DOT);
    plot.setColor(Color.red);
    plot.addPoints(x, g, Plot2.LINE);
    plot.addLabel(0, 0, label);
    if (showApproximation) {
        plot.setColor(Color.blue);
        plot.addPoints(x, f, Plot2.LINE);
    }
    Utils.display(TITLE, plot);
}

Example 2 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class FisherInformationMatrixTest method createFisherInformationMatrix.

private FisherInformationMatrix createFisherInformationMatrix(int n, int k) {
    int maxx = 10;
    int size = maxx * maxx;
    RandomGenerator randomGenerator = new Well19937c(30051977);
    RandomDataGenerator rdg = new RandomDataGenerator(randomGenerator);
    // Use a real Gaussian function here to compute the Fisher information.
    // The matrix may be sensitive to the type of equation used.
    int npeaks = 1;
    while (1 + npeaks * 6 < n) npeaks++;
    Gaussian2DFunction f = GaussianFunctionFactory.create2D(npeaks, maxx, maxx, GaussianFunctionFactory.FIT_ELLIPTICAL, null);
    double[] a = new double[1 + npeaks * 6];
    a[Gaussian2DFunction.BACKGROUND] = rdg.nextUniform(1, 5);
    for (int i = 0, j = 0; i < npeaks; i++, j += 6) {
        a[j + Gaussian2DFunction.SIGNAL] = rdg.nextUniform(100, 300);
        a[j + Gaussian2DFunction.SHAPE] = rdg.nextUniform(-Math.PI, Math.PI);
        // Non-overlapping peaks otherwise the CRLB are poor
        a[j + Gaussian2DFunction.X_POSITION] = rdg.nextUniform(2 + i * 2, 4 + i * 2);
        a[j + Gaussian2DFunction.Y_POSITION] = rdg.nextUniform(2 + i * 2, 4 + i * 2);
        a[j + Gaussian2DFunction.X_SD] = rdg.nextUniform(1.5, 2);
        a[j + Gaussian2DFunction.Y_SD] = rdg.nextUniform(1.5, 2);
    }
    f.initialise(a);
    GradientCalculator c = GradientCalculatorFactory.newCalculator(a.length);
    double[][] I = c.fisherInformationMatrix(size, a, f);
    //System.out.printf("n=%d, k=%d, I=\n", n, k);
    //for (int i = 0; i < I.length; i++)
    //	System.out.println(Arrays.toString(I[i]));
    // Reduce to the desired size
    I = Arrays.copyOf(I, n);
    for (int i = 0; i < n; i++) I[i] = Arrays.copyOf(I[i], n);
    // Zero selected columns
    if (k > 0) {
        int[] zero = new RandomDataGenerator(randomGenerator).nextPermutation(n, k);
        for (int i : zero) {
            for (int j = 0; j < n; j++) {
                I[i][j] = I[j][i] = 0;
            }
        }
    }
    // Create matrix
    return new FisherInformationMatrix(I, 1e-3);
}

Example 3 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian 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;
    }
}

Example 4 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class LVMGradientProcedureTest method gradientProcedureSupportsPrecomputed.

private void gradientProcedureSupportsPrecomputed(final Type type) {
    int iter = 10;
    rdg = new RandomDataGenerator(new Well19937c(30051977));
    ArrayList<double[]> paramsList = new ArrayList<double[]>(iter);
    ArrayList<double[]> yList = new ArrayList<double[]>(iter);
    // 3 peaks
    createData(3, iter, paramsList, yList, true);
    for (int i = 0; i < paramsList.size(); i++) {
        final double[] y = yList.get(i);
        // Add Gaussian read noise so we have negatives
        double min = Maths.min(y);
        for (int j = 0; j < y.length; j++) y[j] = y[i] - min + rdg.nextGaussian(0, Noise);
    }
    // We want to know that:
    // y|peak1+peak2+peak3 == y|peak1+peak2+peak3(precomputed)
    // We want to know when:
    // y|peak1+peak2+peak3 != y-peak3|peak1+peak2
    // i.e. we cannot subtract a precomputed peak from the data, it must be included in the fit
    // E.G. LSQ - subtraction is OK, MLE/WLSQ - subtraction is not allowed
    Gaussian2DFunction f123 = GaussianFunctionFactory.create2D(3, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    Gaussian2DFunction f12 = GaussianFunctionFactory.create2D(2, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    Gaussian2DFunction f3 = GaussianFunctionFactory.create2D(1, blockWidth, blockWidth, GaussianFunctionFactory.FIT_ERF_FREE_CIRCLE, null);
    int nparams = f12.getNumberOfGradients();
    int[] indices = f12.gradientIndices();
    final double[] b = new double[f12.size()];
    double delta = 1e-3;
    DoubleEquality eq = new DoubleEquality(5e-3, 1e-6);
    double[] a1peaks = new double[7];
    final double[] y_b = new double[b.length];
    for (int i = 0; i < paramsList.size(); i++) {
        final double[] y = yList.get(i);
        double[] a3peaks = paramsList.get(i);
        double[] a2peaks = Arrays.copyOf(a3peaks, 1 + 2 * 6);
        double[] a2peaks2 = a2peaks.clone();
        for (int j = 1; j < 7; j++) a1peaks[j] = a3peaks[j + 2 * 6];
        // Evaluate peak 3 to get the background and subtract it from the data to get the new data
        f3.initialise0(a1peaks);
        f3.forEach(new ValueProcedure() {

            int k = 0;

            public void execute(double value) {
                b[k] = value;
                // Remove negatives for MLE
                if (type == Type.MLE) {
                    y[k] = Math.max(0, y[k]);
                    y_b[k] = Math.max(0, y[k] - value);
                } else {
                    y_b[k] = y[k] - value;
                }
                k++;
            }
        });
        // These should be the same
        LVMGradientProcedure p123 = LVMGradientProcedureFactory.create(y, f123, type);
        LVMGradientProcedure p12b3 = LVMGradientProcedureFactory.create(y, PrecomputedGradient1Function.wrapGradient1Function(f12, b), type);
        // This may be different
        LVMGradientProcedure p12m3 = LVMGradientProcedureFactory.create(y_b, f12, type);
        // Check they are the same
        p123.gradient(a3peaks);
        double[][] m123 = p123.getAlphaMatrix();
        p12b3.gradient(a2peaks);
        double s = p12b3.value;
        double[] beta = p12b3.beta.clone();
        double[][] alpha = p12b3.getAlphaMatrix();
        System.out.printf("MLE=%b [%d] p12b3  %f  %f\n", type, i, p123.value, s);
        Assert.assertTrue("p12b3 Not same value @ " + i, eq.almostEqualRelativeOrAbsolute(p123.value, s));
        Assert.assertTrue("p12b3 Not same gradient @ " + i, eq.almostEqualRelativeOrAbsolute(beta, p123.beta));
        for (int j = 0; j < alpha.length; j++) Assert.assertTrue("p12b3 Not same alpha @ " + j, eq.almostEqualRelativeOrAbsolute(alpha[j], m123[j]));
        // Check actual gradients are correct
        for (int j = 0; j < nparams; j++) {
            int k = indices[j];
            double d = Precision.representableDelta(a2peaks[k], (a2peaks[k] == 0) ? 1e-3 : a2peaks[k] * delta);
            a2peaks2[k] = a2peaks[k] + d;
            p12b3.value(a2peaks2);
            double s1 = p12b3.value;
            a2peaks2[k] = a2peaks[k] - d;
            p12b3.value(a2peaks2);
            double s2 = p12b3.value;
            a2peaks2[k] = a2peaks[k];
            // Apply a factor of -2 to compute the actual gradients:
            // See Numerical Recipes in C++, 2nd Ed. Equation 15.5.6 for Nonlinear Models
            beta[j] *= -2;
            double gradient = (s1 - s2) / (2 * d);
            System.out.printf("[%d,%d] %f  (%s %f+/-%f)  %f  ?=  %f  (%f)\n", i, k, s, f12.getName(k), a2peaks[k], d, beta[j], gradient, DoubleEquality.relativeError(gradient, beta[j]));
            Assert.assertTrue("Not same gradient @ " + j, eq.almostEqualRelativeOrAbsolute(beta[j], gradient));
        }
        // Check these may be different
        p12m3.gradient(a2peaks);
        s = p12m3.value;
        beta = p12m3.beta.clone();
        alpha = p12m3.getAlphaMatrix();
        System.out.printf("%s [%d] p12m3  %f  %f\n", type, i, p123.value, s);
        if (type != Type.LSQ) {
            Assert.assertFalse("p12b3 Same value @ " + i, eq.almostEqualRelativeOrAbsolute(p123.value, s));
            Assert.assertFalse("p12b3 Same gradient @ " + i, eq.almostEqualRelativeOrAbsolute(beta, p123.beta));
            for (int j = 0; j < alpha.length; j++) {
                //System.out.printf("%s != %s\n", Arrays.toString(alpha[j]), Arrays.toString(m123[j]));
                Assert.assertFalse("p12b3 Same alpha @ " + j, eq.almostEqualRelativeOrAbsolute(alpha[j], m123[j]));
            }
        } else {
            Assert.assertTrue("p12b3 Not same value @ " + i, eq.almostEqualRelativeOrAbsolute(p123.value, s));
            Assert.assertTrue("p12b3 Not same gradient @ " + i, eq.almostEqualRelativeOrAbsolute(beta, p123.beta));
            for (int j = 0; j < alpha.length; j++) Assert.assertTrue("p12b3 Not same alpha @ " + j, eq.almostEqualRelativeOrAbsolute(alpha[j], m123[j]));
        }
        // Check actual gradients are correct
        for (int j = 0; j < nparams; j++) {
            int k = indices[j];
            double d = Precision.representableDelta(a2peaks[k], (a2peaks[k] == 0) ? 1e-3 : a2peaks[k] * delta);
            a2peaks2[k] = a2peaks[k] + d;
            p12m3.value(a2peaks2);
            double s1 = p12m3.value;
            a2peaks2[k] = a2peaks[k] - d;
            p12m3.value(a2peaks2);
            double s2 = p12m3.value;
            a2peaks2[k] = a2peaks[k];
            // Apply a factor of -2 to compute the actual gradients:
            // See Numerical Recipes in C++, 2nd Ed. Equation 15.5.6 for Nonlinear Models
            beta[j] *= -2;
            double gradient = (s1 - s2) / (2 * d);
            System.out.printf("[%d,%d] %f  (%s %f+/-%f)  %f  ?=  %f  (%f)\n", i, k, s, f12.getName(k), a2peaks[k], d, beta[j], gradient, DoubleEquality.relativeError(gradient, beta[j]));
            Assert.assertTrue("Not same gradient @ " + j, eq.almostEqualRelativeOrAbsolute(beta[j], gradient));
        }
    }
}

Example 5 with Gaussian

use of org.apache.commons.math3.analysis.function.Gaussian in project GDSC-SMLM by aherbert.

the class PCPALMMolecules method addToPlot.

/**
	 * Add the skewed gaussian to the histogram plot
	 * 
	 * @param plot
	 * @param x
	 * @param parameters
	 *            Gaussian parameters
	 * @param alpha
	 * @param shape
	 */
private void addToPlot(Plot2 plot, float[] x, double[] parameters, int shape) {
    SkewNormalFunction sn = new SkewNormalFunction(parameters);
    float[] y = new float[x.length];
    for (int i = 0; i < x.length; i++) y[i] = (float) sn.evaluate(x[i]);
    plot.addPoints(x, y, shape);
}