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

Example 1 with LikelihoodFunction

use of gdsc.smlm.function.LikelihoodFunction 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);
    }
}