Examples with LikelihoodFunction uk.ac.sussex.gdsc.smlm.function.LikelihoodFunction used on opensource projects

Example 1 with LikelihoodFunction

use of uk.ac.sussex.gdsc.smlm.function.LikelihoodFunction in project GDSC-SMLM by aherbert.

the class EmGainAnalysis method plotPmf.

@SuppressWarnings("unused")
private void plotPmf() {
    if (!showPmfDialog()) {
        return;
    }
    final double step = getStepSize(settings.settingPhotons, settings.settingGain, settings.settingNoise);
    final Pdf pdf = pdf(0, step, settings.settingPhotons, settings.settingGain, settings.settingNoise);
    double[] pmf = pdf.probability;
    double yMax = MathUtils.max(pmf);
    // Get the approximation
    LikelihoodFunction fun;
    switch(settings.approximationType) {
        case 3:
            fun = new PoissonFunction(1.0 / settings.settingGain);
            break;
        case 2:
            // Use adaptive normalisation
            fun = PoissonGaussianFunction2.createWithStandardDeviation(1.0 / settings.settingGain, settings.settingNoise);
            break;
        case 1:
            // Create Poisson-Gamma (no Gaussian noise)
            fun = createPoissonGammaGaussianFunction(0);
            break;
        case 0:
        default:
            fun = createPoissonGammaGaussianFunction(settings.settingNoise);
    }
    double expected = settings.settingPhotons;
    if (settings.settingOffset != 0) {
        expected += settings.settingOffset * expected / 100.0;
    }
    // Normalise
    final boolean normalise = false;
    if (normalise) {
        final double sum = MathUtils.sum(pmf);
        for (int i = pmf.length; i-- > 0; ) {
            pmf[i] /= sum;
        }
    }
    // Get CDF
    double sum = 0;
    double sum2 = 0;
    double[] x = pdf.x;
    double[] fvalues = new double[x.length];
    double[] cdf1 = new double[pmf.length];
    double[] cdf2 = new double[pmf.length];
    for (int i = 0; i < cdf1.length; i++) {
        sum += pmf[i] * step;
        cdf1[i] = sum;
        fvalues[i] = fun.likelihood(x[i], expected);
        sum2 += fvalues[i] * step;
        cdf2[i] = sum2;
    }
    // Truncate x for plotting
    int max = 0;
    double plimit = 1 - settings.tail;
    while (sum < plimit && max < pmf.length) {
        sum += pmf[max] * step;
        if (sum > 0.5 && pmf[max] == 0) {
            break;
        }
        max++;
    }
    int min = pmf.length;
    sum = 0;
    plimit = 1 - settings.head;
    while (sum < plimit && min > 0) {
        min--;
        sum += pmf[min] * step;
        if (sum > 0.5 && pmf[min] == 0) {
            break;
        }
    }
    pmf = Arrays.copyOfRange(pmf, min, max);
    x = Arrays.copyOfRange(x, min, max);
    fvalues = Arrays.copyOfRange(fvalues, min, max);
    if (settings.showApproximation) {
        yMax = MathUtils.maxDefault(yMax, fvalues);
    }
    final String label = String.format("Gain=%s, noise=%s, photons=%s", MathUtils.rounded(settings.settingGain), MathUtils.rounded(settings.settingNoise), MathUtils.rounded(settings.settingPhotons));
    final Plot plot = new Plot("PMF", "ADUs", "p");
    plot.setLimits(x[0], x[x.length - 1], 0, yMax);
    plot.setColor(Color.red);
    plot.addPoints(x, pmf, Plot.LINE);
    if (settings.showApproximation) {
        plot.setColor(Color.blue);
        plot.addPoints(x, fvalues, Plot.LINE);
    }
    plot.setColor(Color.magenta);
    plot.drawLine(settings.settingPhotons * settings.settingGain, 0, settings.settingPhotons * settings.settingGain, yMax);
    plot.setColor(Color.black);
    plot.addLabel(0, 0, label);
    final PlotWindow win1 = ImageJUtils.display("PMF", plot);
    // Plot the difference between the actual and approximation
    final double[] delta = new double[fvalues.length];
    for (int i = 0; i < fvalues.length; i++) {
        if (pmf[i] == 0 && fvalues[i] == 0) {
            continue;
        }
        if (settings.relativeDelta) {
            delta[i] = DoubleEquality.relativeError(fvalues[i], pmf[i]) * Math.signum(fvalues[i] - pmf[i]);
        } else {
            delta[i] = fvalues[i] - pmf[i];
        }
    }
    final Plot plot2 = new Plot("PMF delta", "ADUs", (settings.relativeDelta) ? "Relative delta" : "delta");
    final double[] limits = MathUtils.limits(delta);
    plot2.setLimits(x[0], x[x.length - 1], limits[0], limits[1]);
    plot2.setColor(Color.red);
    plot2.addPoints(x, delta, Plot.LINE);
    plot2.setColor(Color.magenta);
    plot2.drawLine(settings.settingPhotons * settings.settingGain, limits[0], settings.settingPhotons * settings.settingGain, limits[1]);
    plot2.setColor(Color.black);
    plot2.addLabel(0, 0, label + ((settings.settingOffset == 0) ? "" : ", expected = " + MathUtils.rounded(expected / settings.settingGain)));
    final WindowOrganiser wo = new WindowOrganiser();
    final PlotWindow win2 = ImageJUtils.display("PMF delta", plot2, wo);
    if (wo.isNotEmpty()) {
        final Point p2 = win1.getLocation();
        p2.y += win1.getHeight();
        win2.setLocation(p2);
    }
    // Plot the CDF of each distribution.
    // Compute the Kolmogorov distance as the supremum (maximum)
    // difference between the two cumulative probability distributions.
    // https://en.wikipedia.org/wiki/Kolmogorov%E2%80%93Smirnov_test
    double kolmogorovDistance = 0;
    double xd = x[0];
    for (int i = 0; i < cdf1.length; i++) {
        final double dist = Math.abs(cdf1[i] - cdf2[i]);
        if (kolmogorovDistance < dist) {
            kolmogorovDistance = dist;
            xd = pdf.x[i];
        }
    }
    cdf1 = Arrays.copyOfRange(cdf1, min, max);
    cdf2 = Arrays.copyOfRange(cdf2, min, max);
    final Plot plot3 = new Plot("CDF", "ADUs", "p");
    yMax = 1.05;
    plot3.setLimits(x[0], x[x.length - 1], 0, yMax);
    plot3.setColor(Color.red);
    plot3.addPoints(x, cdf1, Plot.LINE);
    plot3.setColor(Color.blue);
    plot3.addPoints(x, cdf2, Plot.LINE);
    plot3.setColor(Color.magenta);
    plot3.drawLine(settings.settingPhotons * settings.settingGain, 0, settings.settingPhotons * settings.settingGain, yMax);
    plot3.drawDottedLine(xd, 0, xd, yMax, 2);
    plot3.setColor(Color.black);
    plot3.addLabel(0, 0, label + ", Kolmogorov distance = " + MathUtils.rounded(kolmogorovDistance) + " @ " + xd);
    plot3.addLegend("CDF\nApprox");
    final int size = wo.size();
    final PlotWindow win3 = ImageJUtils.display("CDF", plot3, wo);
    if (size != wo.size()) {
        final Point p2 = win1.getLocation();
        p2.x += win1.getWidth();
        win3.setLocation(p2);
    }
}

Example 2 with LikelihoodFunction

use of uk.ac.sussex.gdsc.smlm.function.LikelihoodFunction in project GDSC-SMLM by aherbert.

the class CameraModelAnalysis method execute.

/**
 * Execute the analysis.
 */
private boolean execute() {
    dirty = false;
    final CameraModelAnalysisSettings settings = this.settings.build();
    if (!(getGain(settings) > 0)) {
        ImageJUtils.log(TITLE + "Error: No total gain");
        return false;
    }
    if (!(settings.getPhotons() > 0)) {
        ImageJUtils.log(TITLE + "Error: No photons");
        return false;
    }
    // Avoid repeating the same analysis
    if (settings.equals(lastSettings)) {
        return true;
    }
    lastSettings = settings;
    final IntHistogram h = getHistogram(settings);
    // Build cumulative distribution
    final double[][] cdf1 = cumulativeHistogram(h);
    final double[] x1 = cdf1[0];
    final double[] y1 = cdf1[1];
    // Interpolate to 300 steps faster evaluation?
    // Get likelihood function
    final LikelihoodFunction f = getLikelihoodFunction(settings);
    // Create likelihood cumulative distribution
    final double[][] cdf2 = cumulativeDistribution(settings, cdf1, f);
    // Compute Komolgorov distance
    final double[] distanceAndValue = getDistance(cdf1, cdf2);
    final double distance = distanceAndValue[0];
    final double value = distanceAndValue[1];
    final double area = distanceAndValue[2];
    final double[] x2 = cdf2[0];
    final double[] y2 = cdf2[1];
    // Fill y1
    int offset = 0;
    while (x2[offset] < x1[0]) {
        offset++;
    }
    final double[] y1b = new double[y2.length];
    System.arraycopy(y1, 0, y1b, offset, y1.length);
    Arrays.fill(y1b, offset + y1.length, y2.length, y1[y1.length - 1]);
    // KolmogorovSmirnovTest
    // n is the number of samples used to build the probability distribution.
    final int n = (int) MathUtils.sum(h.histogramCounts);
    // From KolmogorovSmirnovTest.kolmogorovSmirnovTest(RealDistribution distribution, double[]
    // data, boolean exact):
    // Returns the p-value associated with the null hypothesis that data is a sample from
    // distribution.
    // E.g. If p<0.05 then the null hypothesis is rejected and the data do not match the
    // distribution.
    double pvalue = Double.NaN;
    try {
        pvalue = 1d - kolmogorovSmirnovTest.cdf(distance, n);
    } catch (final MathArithmeticException ex) {
    // Cannot be computed to leave at NaN
    }
    // Plot
    final WindowOrganiser wo = new WindowOrganiser();
    String title = TITLE + " CDF";
    Plot plot = new Plot(title, "Count", "CDF");
    final double max = 1.05 * MathUtils.maxDefault(1, y2);
    plot.setLimits(x2[0], x2[x2.length - 1], 0, max);
    plot.setColor(Color.blue);
    plot.addPoints(x2, y1b, Plot.BAR);
    plot.setColor(Color.red);
    plot.addPoints(x2, y2, Plot.BAR);
    plot.setColor(Color.magenta);
    plot.drawLine(value, 0, value, max);
    plot.setColor(Color.black);
    plot.addLegend("CDF\nModel");
    plot.addLabel(0, 0, String.format("Distance=%s @ %.0f (Mean=%s) : p=%s", MathUtils.rounded(distance), value, MathUtils.rounded(area), MathUtils.rounded(pvalue)));
    ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT, wo);
    // Show the histogram
    title = TITLE + " Histogram";
    plot = new Plot(title, "Count", "Frequency");
    plot.setLimits(x1[0] - 0.5, x1[x1.length - 1] + 1.5, 0, MathUtils.max(h.histogramCounts) * 1.05);
    plot.setColor(Color.blue);
    plot.addPoints(x1, SimpleArrayUtils.toDouble(h.histogramCounts), Plot.BAR);
    plot.setColor(Color.red);
    final double[] x = floatHistogram[0].clone();
    final double[] y = floatHistogram[1].clone();
    final double scale = n / (MathUtils.sum(y) * (x[1] - x[0]));
    for (int i = 0; i < y.length; i++) {
        y[i] *= scale;
    }
    plot.addPoints(x, y, Plot.LINE);
    plot.setColor(Color.black);
    plot.addLegend("Sample\nExpected");
    ImageJUtils.display(title, plot, ImageJUtils.NO_TO_FRONT, wo);
    wo.tile();
    return true;
}