use of org.apache.commons.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);
}
use of org.apache.commons.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);
}
use of org.apache.commons.math3.distribution.PoissonDistribution in project incubator-systemml by apache.
the class PoissonPRNGenerator method setup.
public void setup(double mean, long sd) {
seed = sd;
SynchronizedRandomGenerator srg = new SynchronizedRandomGenerator(new Well1024a());
srg.setSeed(seed);
_pdist = new PoissonDistribution(srg, _mean, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS);
}
use of org.apache.commons.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());
}
use of org.apache.commons.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());
}
Aggregations