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

The following examples show how to use org.apache.commons.math3.distribution.PoissonDistribution#cumulativeProbability() . You can vote up the ones you like or vote down the ones you don't like, and go to the original project or source file by following the links above each example. You may check out the related API usage on the sidebar.
Example 1
Source File: KMap.java    From arx with Apache License 2.0 6 votes vote down vote up
/**
 * 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;
}
 
Example 2
Source File: DriverCatalogFactory.java    From hmftools with GNU General Public License v3.0 5 votes vote down vote up
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()));
}
 
Example 3
Source File: StatisticTest.java    From hmftools with GNU General Public License v3.0 5 votes vote down vote up
@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);

}
 
Example 4
Source File: IndelAnnotator.java    From hmftools with GNU General Public License v3.0 4 votes vote down vote up
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));
            }
        }
    }
}
 
Example 5
Source File: CnJcnCalcs.java    From hmftools with GNU General Public License v3.0 4 votes vote down vote up
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;
}
 
Example 6
Source File: BucketAnalyser.java    From hmftools with GNU General Public License v3.0 4 votes vote down vote up
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));
}
 
Example 7
Source File: CosineSim.java    From hmftools with GNU General Public License v3.0 4 votes vote down vote up
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;
}