Examples with PoissonDistribution uk.ac.sussex.gdsc.smlm.math3.distribution.PoissonDistribution used on opensource projects

Example 36 with PoissonDistribution

use of uk.ac.sussex.gdsc.smlm.math3.distribution.PoissonDistribution in project shifu by ShifuML.

the class AbstractNNWorker method sampleWeights.

protected float sampleWeights(float label) {
    float sampleWeights = 1f;
    // sample negative or kFoldCV, sample rate is 1d
    double sampleRate = (modelConfig.getTrain().getSampleNegOnly() || this.isKFoldCV) ? 1d : modelConfig.getTrain().getBaggingSampleRate();
    int classValue = (int) (label + 0.01f);
    if (!modelConfig.isBaggingWithReplacement()) {
        Random random = null;
        if (this.isStratifiedSampling) {
            random = baggingRandomMap.get(classValue);
            if (random == null) {
                random = DTrainUtils.generateRandomBySampleSeed(modelConfig.getTrain().getBaggingSampleSeed(), CommonConstants.NOT_CONFIGURED_BAGGING_SEED);
                baggingRandomMap.put(classValue, random);
            }
        } else {
            random = baggingRandomMap.get(0);
            if (random == null) {
                random = DTrainUtils.generateRandomBySampleSeed(modelConfig.getTrain().getBaggingSampleSeed(), CommonConstants.NOT_CONFIGURED_BAGGING_SEED);
                baggingRandomMap.put(0, random);
            }
        }
        if (random.nextDouble() <= sampleRate) {
            sampleWeights = 1f;
        } else {
            sampleWeights = 0f;
        }
    } else {
        // replacement
        if (this.isStratifiedSampling) {
            PoissonDistribution rng = this.baggingRngMap.get(classValue);
            if (rng == null) {
                rng = new PoissonDistribution(sampleRate);
                this.baggingRngMap.put(classValue, rng);
            }
            sampleWeights = rng.sample();
        } else {
            PoissonDistribution rng = this.baggingRngMap.get(0);
            if (rng == null) {
                rng = new PoissonDistribution(sampleRate);
                this.baggingRngMap.put(0, rng);
            }
            sampleWeights = rng.sample();
        }
    }
    return sampleWeights;
}

Example 37 with PoissonDistribution

use of uk.ac.sussex.gdsc.smlm.math3.distribution.PoissonDistribution in project GDSC-SMLM by aherbert.

the class PoissonLikelihoodWrapperTest method cumulativeProbabilityIsOneFromLikelihood.

private void cumulativeProbabilityIsOneFromLikelihood(final double mu) {
    // Determine upper limit for a Poisson
    final int limit = new PoissonDistribution(mu).inverseCumulativeProbability(0.999);
    // Expand to allow for the gain
    final int n = (int) Math.ceil(limit / alpha);
    // Evaluate all values from zero to large n
    final double[] k = SimpleArrayUtils.newArray(n, 0, 1.0);
    final double[] a = new double[1];
    final double[] g = new double[1];
    final NonLinearFunction nlf = new DummyNonLinearFunction(mu);
    PoissonLikelihoodWrapper func = new PoissonLikelihoodWrapper(nlf, a, Arrays.copyOf(k, n), n, alpha);
    // Keep evaluating up until no difference
    final double changeTolerance = 1e-6;
    double total = 0;
    double pvalue = 0;
    int index = 0;
    while (index < n) {
        final double nll = func.computeLikelihood(index);
        final double nll2 = func.computeLikelihood(g, index);
        index++;
        Assertions.assertEquals(nll, nll2, 1e-10, "computeLikelihood(i)");
        total += nll;
        final double pp = StdMath.exp(-nll);
        // logger.fine(FunctionUtils.getSupplier("mu=%f, o=%f, i=%d, pp=%f", mu, mu / alpha, i, pp);
        pvalue += pp;
        if (pvalue > 0.5 && pp / pvalue < changeTolerance) {
            break;
        }
    }
    logger.log(TestLogUtils.getRecord(LOG_LEVEL, "mu=%f, limit=%d, p=%f", mu, limit, pvalue));
    Assertions.assertEquals(1, pvalue, 0.02, () -> String.format("mu=%f", mu));
    // Check the function can compute the same total
    func = new PoissonLikelihoodWrapper(nlf, a, k, index, alpha);
    final double sum = func.computeLikelihood();
    final double sum2 = func.computeLikelihood(g);
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
    TestAssertions.assertTest(total, sum, predicate, "computeLikelihood");
    TestAssertions.assertTest(total, sum2, predicate, "computeLikelihood with gradient");
}

Example 38 with PoissonDistribution

use of uk.ac.sussex.gdsc.smlm.math3.distribution.PoissonDistribution in project GDSC-SMLM by aherbert.

the class ScmosLikelihoodWrapperTest method cumulativeProbabilityIsOneWithRealDataForCountAbove8.

@Test
void cumulativeProbabilityIsOneWithRealDataForCountAbove8() {
    for (final double mu : photons) {
        // Determine upper limit for a Poisson
        double max = new PoissonDistribution(mu).inverseCumulativeProbability(P_LIMIT);
        // Determine lower limit
        final double sd = Math.sqrt(mu);
        double min = (int) Math.max(0, mu - 4 * sd);
        // Map to observed values using the gain and offset
        max = max * G + O;
        min = min * G + O;
        cumulativeProbabilityIsOneWithRealData(mu, min, max, mu > 8);
    }
}

Example 39 with PoissonDistribution

use of uk.ac.sussex.gdsc.smlm.math3.distribution.PoissonDistribution in project GDSC-SMLM by aherbert.

the class CameraModelAnalysis method convolveHistogram.

/**
 * Convolve the histogram. The output is a discrete probability distribution.
 *
 * @param settings the settings
 * @return The histogram
 */
private static double[][] convolveHistogram(CameraModelAnalysisSettings settings) {
    // Find the range of the Poisson
    final PoissonDistribution poisson = new PoissonDistribution(settings.getPhotons());
    final int maxn = poisson.inverseCumulativeProbability(UPPER);
    final double gain = getGain(settings);
    final double noise = getReadNoise(settings);
    final boolean debug = false;
    // Build the Probability Mass/Density Function (PDF) of the distribution:
    // either a Poisson (PMF) or Poisson-Gamma (PDF). The PDF is 0 at all
    // values apart from the step interval.
    // Note: The Poisson-Gamma is computed without the Dirac delta contribution
    // at c=0. This allows correct convolution with the Gaussian of the dirac delta
    // and the rest of the Poisson-Gamma (so matching the simulation).
    final TDoubleArrayList list = new TDoubleArrayList();
    double step;
    String name;
    int upsample = 100;
    // Store the Dirac delta value at c=0. This must be convolved separately.
    double dirac = 0;
    // EM-CCD
    if (settings.getMode() == MODE_EM_CCD) {
        name = "Poisson-Gamma";
        final double m = gain;
        final double p = settings.getPhotons();
        dirac = PoissonGammaFunction.dirac(p);
        // Chose whether to compute a discrete PMF or a PDF using the approximation.
        // Note: The delta function at c=0 is from the PMF of the Poisson. So it is
        // a discrete contribution. This is omitted from the PDF and handled in
        // a separate convolution.
        // noise != 0;
        final boolean discrete = false;
        if (discrete) {
            // Note: This is obsolete as the Poisson-Gamma function is continuous.
            // Sampling it at integer intervals is not valid, especially for low gain.
            // The Poisson-Gamma PDF should be integrated to form a discrete PMF.
            step = 1.0;
            double upper;
            if (settings.getPhotons() < 20) {
                upper = maxn;
            } else {
                // Approximate reasonable range of Poisson as a Gaussian
                upper = settings.getPhotons() + 3 * Math.sqrt(settings.getPhotons());
            }
            final GammaDistribution gamma = new GammaDistribution(null, upper, gain, GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
            final int maxc = (int) gamma.inverseCumulativeProbability(0.999);
            final int minn = Math.max(1, poisson.inverseCumulativeProbability(LOWER));
            // See Ulbrich & Isacoff (2007). Nature Methods 4, 319-321, SI equation 3.
            // Note this is not a convolution of a single Gamma distribution since the shape
            // is modified not the count. So it is a convolution of a distribution made with
            // a gamma of fixed count and variable shape.
            // The count=0 is a special case.
            list.add(PoissonGammaFunction.poissonGammaN(0, p, m));
            final long total = (maxn - minn) * (long) maxc;
            if (total < 1000000) {
                // Full computation
                // G(c) = sum n { (1 / n!) p^n e^-p (1 / ((n-1!)m^n)) c^n-1 e^-c/m }
                // Compute as a log
                // - log(n!) + n*log(p)-p -log((n-1)!) - n * log(m) + (n-1) * log(c) -c/m
                // Note: Both methods work
                LogFactorialCache lfc = LOG_FACTORIAL_CACHE.get().get();
                if (lfc == null) {
                    lfc = new LogFactorialCache();
                    LOG_FACTORIAL_CACHE.set(new SoftReference<>(lfc));
                }
                lfc.ensureRange(minn - 1, maxn);
                final double[] f = new double[maxn + 1];
                final double logm = Math.log(m);
                final double logp = Math.log(p);
                for (int n = minn; n <= maxn; n++) {
                    f[n] = -lfc.getLogFactorialUnsafe(n) + n * logp - p - lfc.getLogFactorialUnsafe(n - 1) - n * logm;
                }
                // double total2 = total;
                for (int c = 1; c <= maxc; c++) {
                    double sum = 0;
                    final double c_m = c / m;
                    final double logc = Math.log(c);
                    for (int n = minn; n <= maxn; n++) {
                        sum += StdMath.exp(f[n] + (n - 1) * logc - c_m);
                    }
                    // sum2 += pd[n] * gd[n].density(c);
                    list.add(sum);
                // total += sum;
                // This should match the approximation
                // double approx = PoissonGammaFunction.poissonGamma(c, p, m);
                // total2 += approx;
                // System.out.printf("c=%d sum=%g approx=%g error=%g\n", c, sum2, approx,
                // uk.ac.sussex.gdsc.core.utils.DoubleEquality.relativeError(sum2, approx));
                }
            // System.out.printf("sum=%g approx=%g error=%g\n", total, total2,
            // uk.ac.sussex.gdsc.core.utils.DoubleEquality.relativeError(total, total2));
            } else {
                // Approximate
                for (int c = 1; c <= maxc; c++) {
                    list.add(PoissonGammaFunction.poissonGammaN(c, p, m));
                }
            }
        } else {
            // This integrates the PDF using the approximation and up-samples together.
            // Compute the sampling interval.
            step = 1.0 / upsample;
            // Reset
            upsample = 1;
            // Compute the integral of [-step/2:step/2] for each point.
            // Use trapezoid integration.
            final double step_2 = step / 2;
            double prev = PoissonGammaFunction.poissonGammaN(0, p, m);
            double next = PoissonGammaFunction.poissonGammaN(step_2, p, m);
            list.add((prev + next) * 0.25);
            double max = 0;
            for (int i = 1; ; i++) {
                // Each remaining point is modelling a PMF for the range [-step/2:step/2]
                prev = next;
                next = PoissonGammaFunction.poissonGammaN(i * step + step_2, p, m);
                final double pp = (prev + next) * 0.5;
                if (max < pp) {
                    max = pp;
                }
                if (pp / max < 1e-5) {
                    // Use this if non-zero since it has been calculated
                    if (pp != 0) {
                        list.add(pp);
                    }
                    break;
                }
                list.add(pp);
            }
        }
        // Ensure the combined sum of PDF and Dirac is 1
        final double expected = 1 - dirac;
        // Require an odd number to get an even number (n) of sub-intervals:
        if (list.size() % 2 == 0) {
            list.add(0);
        }
        final double[] g = list.toArray();
        // Number of sub intervals
        final int n = g.length - 1;
        // h = (a-b) / n = sub-interval width
        final double h = 1;
        double sum2 = 0;
        double sum4 = 0;
        for (int j = 1; j <= n / 2 - 1; j++) {
            sum2 += g[2 * j];
        }
        for (int j = 1; j <= n / 2; j++) {
            sum4 += g[2 * j - 1];
        }
        final double sum = (h / 3) * (g[0] + 2 * sum2 + 4 * sum4 + g[n]);
        // Check
        // System.out.printf("Sum=%g Expected=%g\n", sum * step, expected);
        SimpleArrayUtils.multiply(g, expected / sum);
        list.resetQuick();
        list.add(g);
    } else {
        name = "Poisson";
        // Apply fixed gain. Just change the step interval of the PMF.
        step = gain;
        for (int n = 0; n <= maxn; n++) {
            list.add(poisson.probability(n));
        }
        final double p = poisson.probability(list.size());
        if (p != 0) {
            list.add(p);
        }
    }
    // Debug
    if (debug) {
        final String title = name;
        final Plot plot = new Plot(title, "x", "y");
        plot.addPoints(SimpleArrayUtils.newArray(list.size(), 0, step), list.toArray(), Plot.LINE);
        ImageJUtils.display(title, plot);
    }
    double zero = 0;
    double[] pg = list.toArray();
    // Sample Gaussian
    if (noise > 0) {
        step /= upsample;
        pg = list.toArray();
        // Convolve with Gaussian kernel
        final double[] kernel = GaussianKernel.makeGaussianKernel(Math.abs(noise) / step, 6, true);
        if (upsample != 1) {
            // Use scaled convolution. This is faster that zero filling distribution g.
            pg = Convolution.convolve(kernel, pg, upsample);
        } else if (dirac > 0.01) {
            // The Poisson-Gamma may be stepped at low mean causing wrap artifacts in the FFT.
            // This is a problem if most of the probability is in the Dirac.
            // Otherwise it can be ignored and the FFT version is OK.
            pg = Convolution.convolve(kernel, pg);
        } else {
            pg = Convolution.convolveFast(kernel, pg);
        }
        // The convolution will have created a larger array so we must adjust the offset for this
        final int radius = kernel.length / 2;
        zero -= radius * step;
        // Add convolution of the dirac delta function.
        if (dirac != 0) {
            // We only need to convolve the Gaussian at c=0
            for (int i = 0; i < kernel.length; i++) {
                pg[i] += kernel[i] * dirac;
            }
        }
        // Debug
        if (debug) {
            String title = "Gaussian";
            Plot plot = new Plot(title, "x", "y");
            plot.addPoints(SimpleArrayUtils.newArray(kernel.length, -radius * step, step), kernel, Plot.LINE);
            ImageJUtils.display(title, plot);
            title = name + "-Gaussian";
            plot = new Plot(title, "x", "y");
            plot.addPoints(SimpleArrayUtils.newArray(pg.length, zero, step), pg, Plot.LINE);
            ImageJUtils.display(title, plot);
        }
        zero = downSampleCdfToPmf(settings, list, step, zero, pg, 1.0);
        pg = list.toArray();
        zero = (int) Math.floor(zero);
        step = 1.0;
    // No convolution means we have the Poisson PMF/Poisson-Gamma PDF already
    } else if (step != 1) {
        // Sample to 1.0 pixel step interval.
        if (settings.getMode() == MODE_EM_CCD) {
            // Poisson-Gamma PDF
            zero = downSampleCdfToPmf(settings, list, step, zero, pg, 1 - dirac);
            pg = list.toArray();
            zero = (int) Math.floor(zero);
            // Add the dirac delta function.
            if (dirac != 0) {
                // Note: zero is the start of the x-axis. This value should be -1.
                assert (int) zero == -1;
                // Use as an offset to find the actual zero.
                pg[-(int) zero] += dirac;
            }
        } else {
            // Poisson PMF
            // Simple non-interpolated expansion.
            // This should be used when there is no Gaussian convolution.
            final double[] pd = pg;
            list.resetQuick();
            // Account for rounding.
            final Round round = getRound(settings);
            final int maxc = round.round(pd.length * step + 1);
            pg = new double[maxc];
            for (int n = pd.length; n-- > 0; ) {
                pg[round.round(n * step)] += pd[n];
            }
            if (pg[0] != 0) {
                list.add(0);
                list.add(pg);
                pg = list.toArray();
                zero--;
            }
        }
        step = 1.0;
    } else {
        // Add the dirac delta function.
        list.setQuick(0, list.getQuick(0) + dirac);
    }
    return new double[][] { SimpleArrayUtils.newArray(pg.length, zero, step), pg };
}

Example 40 with PoissonDistribution

use of uk.ac.sussex.gdsc.smlm.math3.distribution.PoissonDistribution in project GDSC-SMLM by aherbert.

the class GradientCalculatorSpeedTest method mleGradientCalculatorComputesLikelihood.

@SeededTest
void mleGradientCalculatorComputesLikelihood() {
    // @formatter:off
    final NonLinearFunction func = new NonLinearFunction() {

        double u;

        @Override
        public void initialise(double[] a) {
            u = a[0];
        }

        @Override
        public int[] gradientIndices() {
            return null;
        }

        @Override
        public double eval(int x, double[] dyda) {
            return 0;
        }

        @Override
        public double eval(int x) {
            return u;
        }

        @Override
        public double evalw(int x, double[] dyda, double[] w) {
            return 0;
        }

        @Override
        public double evalw(int x, double[] w) {
            return 0;
        }

        @Override
        public boolean canComputeWeights() {
            return false;
        }

        @Override
        public int getNumberOfGradients() {
            return 0;
        }
    };
    // @formatter:on
    final DoubleDoubleBiPredicate predicate = TestHelper.doublesAreClose(1e-10, 0);
    final int[] xx = SimpleArrayUtils.natural(100);
    final double[] xxx = SimpleArrayUtils.newArray(100, 0, 1.0);
    for (final double u : new double[] { 0.79, 2.5, 5.32 }) {
        double ll = 0;
        double oll = 0;
        final PoissonDistribution pd = new PoissonDistribution(u);
        // The logLikelihood function for the entire array of observations is then asserted.
        for (final int x : xx) {
            double obs = PoissonCalculator.likelihood(u, x);
            double exp = pd.probability(x);
            TestAssertions.assertTest(exp, obs, predicate, "likelihood");
            obs = PoissonCalculator.logLikelihood(u, x);
            exp = pd.logProbability(x);
            TestAssertions.assertTest(exp, obs, predicate, "log likelihood");
            oll += obs;
            ll += exp;
        }
        final MleGradientCalculator gc = new MleGradientCalculator(1);
        final double o = gc.logLikelihood(xxx, new double[] { u }, func);
        Assertions.assertEquals(oll, o, "sum log likelihood should exactly match the PoissonCalculator");
        TestAssertions.assertTest(ll, o, predicate, "sum log likelihood");
    }
}