Java Code Examples for org.apache.commons.math3.distribution.PoissonDistribution#cumulativeProbability()

/**
 * Calculates k, based on Poisson distribution.
 * @param lambda
 * @return
 */
private int calculateKPoisson(double lambda) {
    
    final double threshold = 1d - this.significanceLevel;
    final PoissonDistribution distribution = new PoissonDistribution(lambda);
    int counter = 0;
    double value = 0;
    while (value < threshold) {
        // value += (Math.pow(lambda, counter) * Math.exp(-lambda)) / ArithmeticUtils.factorial(counter);
        value = distribution.cumulativeProbability(counter);
        counter++;
        // Break if the estimated k is equal or greater than the given k, as this makes no sense.
        if (counter >= this.k) {
            // We are 100% sure that the dataset fulfills k-map
            value = 1d;
            break;
        }
    }
    this.type1Error = 1d - value;
    return counter + 1;
}
 

private static double probabilityDriverVariantSameImpact(int count, long sampleSNVCount,
        @NotNull final DndsDriverImpactLikelihood likelihood) {
    double lambda = sampleSNVCount * likelihood.pVariantNonDriverFactor();
    if (Doubles.isZero(lambda)) {
        return likelihood.dndsLikelihood();
    }

    PoissonDistribution poissonDistribution = new PoissonDistribution(lambda);

    double pVariantNonDriver = 1 - poissonDistribution.cumulativeProbability(count);
    return likelihood.pDriver() / (likelihood.pDriver() + pVariantNonDriver * (1 - likelihood.pDriver()));
}
 

@Test
public void testExternalFunctions()
{
    // chiSquaredTests
    int degreesOfFreedom = 95;
    ChiSquaredDistribution chiSquDist = new ChiSquaredDistribution(degreesOfFreedom);

    double result = chiSquDist.cumulativeProbability(135);

    PoissonDistribution poisson = new PoissonDistribution(100);

    double prob = poisson.cumulativeProbability(99);
    prob = poisson.cumulativeProbability(110);

}
 

public void annotateCluster(final SvCluster cluster)
{
    if(cluster.getResolvedType() != COMPLEX && cluster.getResolvedType() != SIMPLE_GRP && cluster.getResolvedType() != RECIP_INV)
        return;

    int totalIndels = 0;

    for(SvLinkedPair pair : cluster.getLinkedPairs())
    {
        int lowerPos = pair.getBreakend(true).position();
        int upperPos = pair.getBreakend(false).position();
        pair.setIndelCount(findIndels(pair.chromosome(), lowerPos, upperPos).size());

        if(pair.getIndelCount() > 0)
        {
            totalIndels += pair.getIndelCount();
            LNX_LOGGER.debug("cluster({}) pair({} len={}) indelCount({})",
                    cluster.id(), pair, pair.length(), pair.getIndelCount());
        }
    }

    ClusterMetrics metrics = cluster.getMetrics();
    metrics.IndelCount = totalIndels;
    metrics.IndelProbability = 1; // default indicates nothing unexpected

    if(!cluster.hasVariedJcn() && metrics.TotalRange > 0 && metrics.TotalDeleted < metrics.TotalRange)
    {
        long clusterRange = metrics.TotalRange - metrics.TotalDeleted;
        double expectedIndelCount = min(clusterRange / GENOME_LENGTH, 1) * (mIndelCount - mMsiIndelCount);

        if (totalIndels > expectedIndelCount)
        {
            PoissonDistribution poisson = new PoissonDistribution(expectedIndelCount);
            double indelProbability = 1 - poisson.cumulativeProbability(totalIndels - 1);
            metrics.IndelProbability = indelProbability;

            if(indelProbability < 0.01)
            {
                LNX_LOGGER.debug(String.format("cluster(%d) indelCount(%d) range(%d) expected(%.1f) probability(%.9f)",
                        cluster.id(), totalIndels, clusterRange, expectedIndelCount, indelProbability));
            }
        }
    }
}
 

private static int calcPoissonObservedGivenProb(int expectedVal, double requiredProb)
{
    if(expectedVal <= 0)
        return 0;

    PoissonDistribution poisson = new PoissonDistribution(expectedVal);

    int maxIterations = 20;
    int iterations = 0;

    double refCount = 25;
    double refRangePerc = 0.44;
    double rangePerc = refRangePerc / sqrt(expectedVal/refCount);
    double range = rangePerc * expectedVal;

    int testValueUpper;
    int testValueLower;

    if(requiredProb > 0.5)
    {
        testValueUpper = (int) round(expectedVal + range * 1.5);
        testValueLower = (int) round(max(expectedVal + range * 0.2, 0));
    }
    else
    {
        testValueUpper = (int) round(expectedVal - range * 0.2);
        testValueLower = (int) round(max(expectedVal - range * 1.2, 0));
    }

    int testValue = (int)((testValueLower + testValueUpper) * 0.5);

    double currentProb = poisson.cumulativeProbability(testValue);
    double probDiff = 0;

    while(iterations < maxIterations)
    {
        probDiff = abs(requiredProb - currentProb) / requiredProb;

        if(probDiff < 0.001)
            break;

        // if prob is too high, need to lower the test value
        if(currentProb > requiredProb)
        {
            if(testValue <= testValueLower + 1)
                break;

            testValueUpper = testValue;
            testValue = (int)round((testValue + testValueLower) * 0.5);
        }
        else
        {
            if(testValue >= testValueUpper - 1)
                break;

            testValueLower = testValue;
            testValue = (int)round((testValue + testValueUpper) * 0.5);
        }

        currentProb = poisson.cumulativeProbability(testValue);
        ++iterations;
    }

    if(iterations >= maxIterations)
    {
        LNX_LOGGER.warn(String.format("max iterations reached: value(%d) test(%d) prob(%.4f diff=%.4f)",
                expectedVal, testValue, currentProb, probDiff));
    }

    return testValue;
}
 

private void splitSampleCounts()
{
    // split each sample's bucket counts into background and elevated
    // and additionally calculate a probability for each bucket that it is elevated above the background
    LOGGER.debug("splitting sample counts");

    // work out bucket median values (literally 50th percentile values)
    mBucketProbs = new SigMatrix(mBucketCount, mSampleCount);
    mElevatedCounts = new SigMatrix(mBucketCount, mSampleCount);

    double[][] probData = mBucketProbs.getData();
    double[][] bgData = mBackgroundSigDiscovery.getBackgroundCounts().getData();
    double[][] elevData = mElevatedCounts.getData();
    final double[][] scData = mSampleCounts.getData();

    // record the frequency of each calculated probability - purely for logging purposes
    int probExpMax = 15;
    int[] probFrequencies = new int[probExpMax+1];
    double minProb = pow(10, -probExpMax);
    int zeroProbIndex = probFrequencies.length-1;

    int gridSize = mSampleCount * mBucketCount;

    for (int i = 0; i < mBucketCount; ++i)
    {
        for(int j = 0; j < mSampleCount; ++j)
        {
            int sbCount = (int)scData[i][j];

            int backgroundCount = (int)bgData[i][j];

            if(sbCount < backgroundCount) // must be capped by the actual count
                backgroundCount = sbCount;

            int elevatedCount = sbCount - backgroundCount;

            // compute a range for Poisson noise around this elevated count
            if (elevatedCount > 0)
            {
                elevData[i][j] = elevatedCount;
            }

            double prob = 1;

            if (backgroundCount == 0)
            {
                prob = 0;
            }
            else if (elevatedCount > 0)
            {
                PoissonDistribution poisDist = new PoissonDistribution(backgroundCount);
                prob = 1 - poisDist.cumulativeProbability(sbCount - 1);
                prob = min(prob * gridSize, 1); // apply false discovery rate being # tests
            }

            probData[i][j] = prob;

            if (prob > minProb && prob < 0.1)
            {
                int baseProb = -(int) round(log10(prob));

                if (baseProb >= 0 && baseProb <= probExpMax)
                    probFrequencies[baseProb] += 1;
            }
            else if (prob < minProb)
            {
                // allocate to the last slot
                probFrequencies[zeroProbIndex] += 1;
            }
        }
    }

    mElevatedCounts.cacheTranspose();

    mBackgroundCount = mBackgroundSigDiscovery.getBackgroundCounts().sum();
    mElevatedCount = mElevatedCounts.sum();

    LOGGER.debug(String.format("total counts: background(%s perc=%.3f) elevated(%s perc=%.3f) of total(%s)",
            sizeToStr(mBackgroundCount), mBackgroundCount/mTotalCount,
            sizeToStr(mElevatedCount), mElevatedCount/mTotalCount, sizeToStr(mTotalCount)));

    for (int i = 0; i < probFrequencies.length-1; ++i)
    {
        LOGGER.debug(String.format("probability(1e-%d) freq(%d) percOfTotal(%.4f)",
                i, probFrequencies[i], probFrequencies[i]/(double)gridSize));
    }

    LOGGER.debug(String.format("probability(zero) freq(%d) percOfTotal(%.4f)",
            probFrequencies[zeroProbIndex], probFrequencies[zeroProbIndex]/(double)gridSize));
}
 

public static int calcPoissonRangeGivenProb(int value, double requiredProb)
{
    if(value <= 10)
        return 10;

    PoissonDistribution poisson = new PoissonDistribution(value);

    int maxIterations = 10;
    int iterations = 0;

    double initRange = 3.7 / sqrt(value); // works for requiredProb = 1e-4
    int testValue = (int)max(round(value * (1 - initRange)), 0);
    int testValueUpper = (int)max(round(value * (1 - initRange*0.5)), 0);
    int testValueLower = (int)max(round(value * (1 - initRange*2)), 0);

    double currentProb = poisson.cumulativeProbability(testValue);
    double probDiff = 0;

    while(iterations < maxIterations)
    {
        probDiff = abs(requiredProb - currentProb) / requiredProb;

        if(probDiff < 0.1)
            break;

        // if prob is too high, need to lower the test value
        if(currentProb > requiredProb)
        {
            if(testValue <= testValueLower + 1)
                break;

            testValueUpper = testValue;
            testValue = (int)round((testValue + testValueLower) * 0.5);
        }
        else
        {
            if(testValue >= testValueUpper - 1)
                break;

            testValueLower = testValue;
            testValue = (int)round((testValue + testValueUpper) * 0.5);
        }

        currentProb = poisson.cumulativeProbability(testValue);
        ++iterations;
    }

    if(iterations >= maxIterations)
    {
        LOGGER.warn(String.format("max iterations reached: value(%d) test(%d) prob(%.4f diff=%.4f)", value, testValue, currentProb, probDiff));
    }

    return value - testValue;
}