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

Example 21 with PoissonDistribution

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

the class SCMOSLikelihoodWrapperTest method instanceLikelihoodMatches.

private void instanceLikelihoodMatches(final double mu, boolean test) {
    // Determine upper limit for a Poisson
    int limit = new PoissonDistribution(mu).inverseCumulativeProbability(P_LIMIT);
    // Map to observed values using the gain and offset
    double max = limit * G;
    double step = 0.1;
    int n = (int) Math.ceil(max / step);
    // Evaluate all values from (zero+offset) to large n
    double[] k = Utils.newArray(n, O, step);
    double[] a = new double[0];
    double[] gradient = new double[0];
    float[] var = newArray(n, VAR);
    float[] g = newArray(n, G);
    float[] o = newArray(n, O);
    NonLinearFunction nlf = new NonLinearFunction() {

        public void initialise(double[] a) {
        }

        public int[] gradientIndices() {
            return new int[0];
        }

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

        public double eval(int x) {
            return mu;
        }

        public double eval(int x, double[] dyda) {
            return mu;
        }

        public boolean canComputeWeights() {
            return false;
        }

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

        public int getNumberOfGradients() {
            return 0;
        }
    };
    SCMOSLikelihoodWrapper f = new SCMOSLikelihoodWrapper(nlf, a, k, n, var, g, o);
    double total = 0, p = 0;
    double maxp = 0;
    int maxi = 0;
    for (int i = 0; i < n; i++) {
        double nll = f.computeLikelihood(i);
        double nll2 = f.computeLikelihood(gradient, i);
        double nll3 = SCMOSLikelihoodWrapper.negativeLogLikelihood(mu, var[i], g[i], o[i], k[i]);
        total += nll;
        Assert.assertEquals("computeLikelihood @" + i, nll3, nll, nll * 1e-10);
        Assert.assertEquals("computeLikelihood+gradient @" + i, nll3, nll2, nll * 1e-10);
        double pp = FastMath.exp(-nll);
        if (maxp < pp) {
            maxp = pp;
            maxi = i;
        //System.out.printf("mu=%f, e=%f, k=%f, pp=%f\n", mu, mu * G + O, k[i], pp);
        }
        p += pp * step;
    }
    // Expected max of the distribution is the mode of the Poisson distribution.
    // This has two modes for integer input counts. We take the mean of those.
    // https://en.wikipedia.org/wiki/Poisson_distribution
    // Note that the shift of VAR/(G*G) is a constant applied to both the expected and
    // observed values and consequently cancels when predicting the max, i.e. we add
    // a constant count to the observed values and shift the distribution by the same
    // constant. We can thus compute the mode for the unshifted distribution.
    double lambda = mu;
    double mode1 = Math.floor(lambda);
    double mode2 = Math.ceil(lambda) - 1;
    // Scale to observed values
    double kmax = ((mode1 + mode2) * 0.5) * G + O;
    //System.out.printf("mu=%f, p=%f, maxp=%f @ %f  (expected=%f  %f)\n", mu, p, maxp, k[maxi], kmax, kmax - k[maxi]);
    Assert.assertEquals("k-max", kmax, k[maxi], kmax * 1e-3);
    if (test) {
        Assert.assertEquals(String.format("mu=%f", mu), P_LIMIT, p, 0.02);
    }
    // Check the function can compute the same total
    double delta = total * 1e-10;
    double sum, sum2;
    sum = f.computeLikelihood();
    sum2 = f.computeLikelihood(gradient);
    Assert.assertEquals("computeLikelihood", total, sum, delta);
    Assert.assertEquals("computeLikelihood with gradient", total, sum2, delta);
    // Check the function can compute the same total after duplication
    f = f.build(nlf, a);
    sum = f.computeLikelihood();
    sum2 = f.computeLikelihood(gradient);
    Assert.assertEquals("computeLikelihood", total, sum, delta);
    Assert.assertEquals("computeLikelihood with gradient", total, sum2, delta);
}

Example 22 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
    int limit = new PoissonDistribution(mu).inverseCumulativeProbability(0.999);
    // Expand to allow for the gain
    int n = (int) Math.ceil(limit / alpha);
    // Evaluate all values from zero to large n
    double[] k = Utils.newArray(n, 0, 1.0);
    double[] a = new double[0];
    double[] g = new double[0];
    NonLinearFunction nlf = new NonLinearFunction() {

        public void initialise(double[] a) {
        }

        public int[] gradientIndices() {
            return new int[0];
        }

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

        public double eval(int x) {
            return mu / alpha;
        }

        public double eval(int x, double[] dyda) {
            return mu / alpha;
        }

        public boolean canComputeWeights() {
            return false;
        }

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

        public int getNumberOfGradients() {
            return 0;
        }
    };
    PoissonLikelihoodWrapper f = 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 p = 0;
    int i = 0;
    while (i < n) {
        double nll = f.computeLikelihood(i);
        double nll2 = f.computeLikelihood(g, i);
        i++;
        Assert.assertEquals("computeLikelihood(i)", nll, nll2, 1e-10);
        total += nll;
        double pp = FastMath.exp(-nll);
        //System.out.printf("mu=%f, o=%f, i=%d, pp=%f\n", mu, mu / alpha, i, pp);
        p += pp;
        if (p > 0.5 && pp / p < changeTolerance)
            break;
    }
    System.out.printf("mu=%f, limit=%d, p=%f\n", mu, limit, p);
    Assert.assertEquals(String.format("mu=%f", mu), 1, p, 0.02);
    // Check the function can compute the same total
    f = new PoissonLikelihoodWrapper(nlf, a, k, i, alpha);
    double sum = f.computeLikelihood();
    double sum2 = f.computeLikelihood(g);
    double delta = total * 1e-10;
    Assert.assertEquals("computeLikelihood", total, sum, delta);
    Assert.assertEquals("computeLikelihood with gradient", total, sum2, delta);
}

Example 23 with PoissonDistribution

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

the class EMGainAnalysis method simulateFromPoissonGammaGaussian.

/**
	 * Randomly generate a histogram from poisson-gamma-gaussian samples
	 * 
	 * @return The histogram
	 */
private int[] simulateFromPoissonGammaGaussian() {
    // Randomly sample
    RandomGenerator random = new Well44497b(System.currentTimeMillis() + System.identityHashCode(this));
    PoissonDistribution poisson = new PoissonDistribution(random, _photons, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
    CustomGammaDistribution gamma = new CustomGammaDistribution(random, _photons, _gain, GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
    final int steps = simulationSize;
    int[] sample = new int[steps];
    for (int n = 0; n < steps; n++) {
        if (n % 64 == 0)
            IJ.showProgress(n, steps);
        // Poisson
        double d = poisson.sample();
        // Gamma
        if (d > 0) {
            gamma.setShapeUnsafe(d);
            d = gamma.sample();
        }
        // Gaussian
        d += _noise * random.nextGaussian();
        // Convert the sample to a count 
        sample[n] = (int) Math.round(d + _bias);
    }
    int max = Maths.max(sample);
    int[] h = new int[max + 1];
    for (int s : sample) h[s]++;
    return h;
}

Example 24 with PoissonDistribution

use of uk.ac.sussex.gdsc.smlm.math3.distribution.PoissonDistribution in project gatk by broadinstitute.

the class CoverageModelCopyRatioEmissionProbabilityCalculator method logLikelihoodPoisson.

/**
     * Calculate emission log probability directly using Poisson distribution
     *
     * TODO github/gatk-protected issue #854
     *
     * @implNote This is a naive implementation where the variance of log bias ($\Psi$) is
     * ignored altogether. This routine must be improved by performing a few-point numerical
     * integration of:
     *
     *      \int_{-\infty}^{+\infty} db Poisson(n | \lambda = d*c*exp(b) + eps*d)
     *                                                           * Normal(b | \mu, \Psi)
     *
     * @param emissionData an instance of {@link CoverageModelCopyRatioEmissionData}
     * @param copyRatio copy ratio on which the emission probability is conditioned on
     * @return a double value
     */
private double logLikelihoodPoisson(@Nonnull CoverageModelCopyRatioEmissionData emissionData, double copyRatio) {
    final double multBias = FastMath.exp(emissionData.getMu());
    final double readDepth = emissionData.getCopyRatioCallingMetadata().getSampleCoverageDepth();
    final double err = emissionData.getMappingErrorProbability();
    final double poissonMean = readDepth * ((1 - err) * copyRatio * multBias + err);
    return new PoissonDistribution(null, poissonMean, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS).logProbability(emissionData.getReadCount());
}

Example 25 with PoissonDistribution

use of uk.ac.sussex.gdsc.smlm.math3.distribution.PoissonDistribution in project gatk by broadinstitute.

the class TargetCoverageSexGenotypeCalculator method calculateSexGenotypeData.

/**
     * Calculates the likelihood of different sex genotypes for a given sample index
     *
     * @param sampleIndex sample index
     * @return an instance of {@link SexGenotypeData}
     */
private SexGenotypeData calculateSexGenotypeData(final int sampleIndex) {
    final int[] allosomalReadCounts = Arrays.stream(processedReadCounts.getColumnOnSpecifiedTargets(sampleIndex, allosomalTargetList, true)).mapToInt(n -> (int) n).toArray();
    final double readDepth = getSampleReadDepthFromAutosomalTargets(sampleIndex);
    logger.debug("Read depth for " + processedReadCounts.columnNames().get(sampleIndex) + ": " + readDepth);
    final List<Double> logLikelihoods = new ArrayList<>();
    final List<String> sexGenotypesList = new ArrayList<>();
    final int numAllosomalTargets = allosomalTargetList.size();
    for (final String genotypeName : sexGenotypeIdentifiers) {
        /* calculate log likelihood */
        final int[] currentAllosomalTargetPloidies = allosomalTargetPloidies.get(genotypeName);
        final double[] poissonParameters = IntStream.range(0, numAllosomalTargets).mapToDouble(i -> readDepth * (currentAllosomalTargetPloidies[i] > 0 ? currentAllosomalTargetPloidies[i] : baselineMappingErrorProbability)).toArray();
        final double currentLogLikelihood = IntStream.range(0, numAllosomalTargets).mapToDouble(i -> {
            final PoissonDistribution pois = new PoissonDistribution(poissonParameters[i]);
            return pois.logProbability(allosomalReadCounts[i]);
        }).sum();
        sexGenotypesList.add(genotypeName);
        logLikelihoods.add(currentLogLikelihood);
    }
    /* infer the most likely sex genotype */
    final Integer[] indices = new Integer[sexGenotypesList.size()];
    IntStream.range(0, sexGenotypesList.size()).forEach(i -> indices[i] = i);
    Arrays.sort(indices, (li, ri) -> Double.compare(logLikelihoods.get(ri), logLikelihoods.get(li)));
    return new SexGenotypeData(processedReadCounts.columnNames().get(sampleIndex), sexGenotypesList.get(indices[0]), sexGenotypesList, logLikelihoods.stream().mapToDouble(d -> d).toArray());
}