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 |
/** * 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 |
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 |
@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 |
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 |
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 |
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 |
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; }