Examples with PoissonDistributionImpl org.apache.commons.math.distribution.PoissonDistributionImpl used on opensource projects

Example 1 with PoissonDistributionImpl

use of org.apache.commons.math.distribution.PoissonDistributionImpl in project bacter by tgvaughan.

the class ACGCoalescent method calculateLogP.

@Override
public double calculateLogP() {
    // Check whether conversion count exceeds bounds.
    if (acg.getTotalConvCount() < lowerCCBoundInput.get() || acg.getTotalConvCount() > upperCCBoundInput.get())
        return Double.NEGATIVE_INFINITY;
    logP = calculateClonalFrameLogP();
    double poissonMean = rhoInput.get().getValue() * acg.getClonalFrameLength() * (acg.getTotalSequenceLength() + acg.getLoci().size() * (deltaInput.get().getValue() - 1.0));
    // Probability of conversion count:
    if (poissonMean > 0.0) {
        logP += -poissonMean + acg.getTotalConvCount() * Math.log(poissonMean);
    // - GammaFunction.lnGamma(acg.getConvCount()+1);
    } else {
        if (acg.getTotalConvCount() > 0)
            logP = Double.NEGATIVE_INFINITY;
    }
    for (Locus locus : acg.getLoci()) for (Conversion conv : acg.getConversions(locus)) logP += calculateConversionLogP(conv);
    if (lowerCCBoundInput.get() > 0 || upperCCBoundInput.get() < Integer.MAX_VALUE) {
        try {
            logP -= new PoissonDistributionImpl(poissonMean).cumulativeProbability(lowerCCBoundInput.get(), upperCCBoundInput.get());
        } catch (MathException e) {
            throw new RuntimeException("Error computing modification to ARG " + "prior density required by conversion number constraint.");
        }
    }
    return logP;
}

Example 2 with PoissonDistributionImpl

use of org.apache.commons.math.distribution.PoissonDistributionImpl in project beast-mcmc by beast-dev.

the class PoissonDistribution method logPdf.

public static double logPdf(double x, double mean) {
    PoissonDistributionImpl dist = new PoissonDistributionImpl(mean);
    double pdf = dist.probability(x);
    if (pdf == 0 || Double.isNaN(pdf)) {
        // bad estimate
        return x * Math.log(mean) - Poisson.gammln(x + 1) - mean;
    }
    return Math.log(pdf);
}