Examples with PoissonGammaGaussianFunction gdsc.smlm.function.PoissonGammaGaussianFunction used on opensource projects

Example 1 with PoissonGammaGaussianFunction

use of gdsc.smlm.function.PoissonGammaGaussianFunction 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 PoissonGammaGaussianFunction

use of gdsc.smlm.function.PoissonGammaGaussianFunction in project GDSC-SMLM by aherbert.

the class EMGainAnalysis method plotPMF.

private void plotPMF() {
    if (!showPMFDialog())
        return;
    final int gaussWidth = 5;
    int dummyBias = (int) Math.max(500, gaussWidth * _noise + 1);
    double[] pmf = pdf(0, _photons, _gain, _noise, dummyBias);
    double[] x = Utils.newArray(pmf.length, 0, 1.0);
    double yMax = Maths.max(pmf);
    // Truncate x
    int max = 0;
    double sum = 0;
    double p = 1 - tail;
    while (sum < p && max < pmf.length) {
        sum += pmf[max];
        if (sum > 0.5 && pmf[max] == 0)
            break;
        max++;
    }
    int min = pmf.length;
    sum = 0;
    p = 1 - head;
    while (sum < p && min > 0) {
        min--;
        sum += pmf[min];
        if (sum > 0.5 && pmf[min] == 0)
            break;
    }
    //int min = (int) (dummyBias - gaussWidth * _noise);
    pmf = Arrays.copyOfRange(pmf, min, max);
    x = Arrays.copyOfRange(x, min, max);
    // Get the approximation
    double[] f = new double[x.length];
    LikelihoodFunction fun;
    double myNoise = _noise;
    switch(approximation) {
        case 3:
            fun = new PoissonFunction(1.0 / _gain, true);
            break;
        case 2:
            // The mean does not matter so just use zero
            fun = PoissonGaussianFunction.createWithStandardDeviation(1.0 / _gain, 0, _noise);
            break;
        case 1:
            myNoise = 0;
        case 0:
        default:
            PoissonGammaGaussianFunction myFun = new PoissonGammaGaussianFunction(1.0 / _gain, myNoise);
            myFun.setMinimumProbability(0);
            fun = myFun;
    }
    double expected = _photons;
    if (offset != 0)
        expected += offset * expected / 100.0;
    expected *= _gain;
    //double sum2 = 0;
    for (int i = 0; i < f.length; i++) {
        // Adjust the x-values to remove the dummy bias
        x[i] -= dummyBias;
        f[i] = fun.likelihood(x[i], expected);
    //sum += pmf[i];
    //sum2 += f[i];
    }
    //System.out.printf("Approximation sum = %f : %f\n", sum ,sum2);
    if (showApproximation)
        yMax = Maths.maxDefault(yMax, f);
    String label = String.format("Gain=%s, noise=%s, photons=%s", Utils.rounded(_gain), Utils.rounded(_noise), Utils.rounded(_photons));
    Plot2 plot = new Plot2("PMF", "ADUs", "p");
    plot.setLimits(x[0], x[x.length - 1], 0, yMax);
    plot.setColor(Color.red);
    plot.addPoints(x, pmf, Plot2.LINE);
    if (showApproximation) {
        plot.setColor(Color.blue);
        plot.addPoints(x, f, Plot2.LINE);
    }
    plot.setColor(Color.magenta);
    plot.drawLine(_photons * _gain, 0, _photons * _gain, yMax);
    plot.setColor(Color.black);
    plot.addLabel(0, 0, label);
    PlotWindow win1 = Utils.display("PMF", plot);
    // Plot the difference between the actual and approximation
    double[] delta = new double[f.length];
    for (int i = 0; i < f.length; i++) {
        if (pmf[i] == 0 && f[i] == 0)
            continue;
        if (relativeDelta)
            delta[i] = DoubleEquality.relativeError(f[i], pmf[i]) * Math.signum(f[i] - pmf[i]);
        else
            delta[i] = f[i] - pmf[i];
    }
    Plot2 plot2 = new Plot2("PMF delta", "ADUs", (relativeDelta) ? "Relative delta" : "delta");
    double[] limits = Maths.limits(delta);
    plot2.setLimits(x[0], x[x.length - 1], limits[0], limits[1]);
    plot2.setColor(Color.red);
    plot2.addPoints(x, delta, Plot2.LINE);
    plot2.setColor(Color.magenta);
    plot2.drawLine(_photons * _gain, limits[0], _photons * _gain, limits[1]);
    plot2.setColor(Color.black);
    plot2.addLabel(0, 0, label + ((offset == 0) ? "" : ", expected = " + Utils.rounded(expected / _gain)));
    PlotWindow win2 = Utils.display("PMF delta", plot2);
    if (Utils.isNewWindow()) {
        Point p2 = win2.getLocation();
        p2.y += win1.getHeight();
        win2.setLocation(p2);
    }
}