org.apache.commons.math3.distribution.BinomialDistribution Java Examples

@Test
public void testNextBinomial() {
    BinomialDistributionTest testInstance = new BinomialDistributionTest();
    int[] densityPoints = testInstance.makeDensityTestPoints();
    double[] densityValues = testInstance.makeDensityTestValues();
    int sampleSize = 1000;
    int length = TestUtils.eliminateZeroMassPoints(densityPoints, densityValues);
    BinomialDistribution distribution = (BinomialDistribution) testInstance.makeDistribution();
    double[] expectedCounts = new double[length];
    long[] observedCounts = new long[length];
    for (int i = 0; i < length; i++) {
        expectedCounts[i] = sampleSize * densityValues[i];
    }
    randomData.reSeed(1000);
    for (int i = 0; i < sampleSize; i++) {
      int value = randomData.nextBinomial(distribution.getNumberOfTrials(),
              distribution.getProbabilityOfSuccess());
      for (int j = 0; j < length; j++) {
          if (value == densityPoints[j]) {
              observedCounts[j]++;
          }
      }
    }
    TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, observedCounts, .001);
}
 

@Test
public void testNextBinomial() {
    BinomialDistributionTest testInstance = new BinomialDistributionTest();
    int[] densityPoints = testInstance.makeDensityTestPoints();
    double[] densityValues = testInstance.makeDensityTestValues();
    int sampleSize = 1000;
    int length = TestUtils.eliminateZeroMassPoints(densityPoints, densityValues);
    BinomialDistribution distribution = (BinomialDistribution) testInstance.makeDistribution();
    double[] expectedCounts = new double[length];
    long[] observedCounts = new long[length];
    for (int i = 0; i < length; i++) {
        expectedCounts[i] = sampleSize * densityValues[i];
    }
    randomData.reSeed(1000);
    for (int i = 0; i < sampleSize; i++) {
      int value = randomData.nextBinomial(distribution.getNumberOfTrials(),
              distribution.getProbabilityOfSuccess());
      for (int j = 0; j < length; j++) {
          if (value == densityPoints[j]) {
              observedCounts[j]++;
          }
      }
    }
    TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, observedCounts, .001);
}
 

@Test
public void testNextBinomial() throws Exception {
    BinomialDistributionTest testInstance = new BinomialDistributionTest();
    int[] densityPoints = testInstance.makeDensityTestPoints();
    double[] densityValues = testInstance.makeDensityTestValues();
    int sampleSize = 1000;
    int length = TestUtils.eliminateZeroMassPoints(densityPoints, densityValues);
    BinomialDistribution distribution = (BinomialDistribution) testInstance.makeDistribution();
    double[] expectedCounts = new double[length];
    long[] observedCounts = new long[length];
    for (int i = 0; i < length; i++) {
        expectedCounts[i] = sampleSize * densityValues[i];
    }
    randomData.reSeed(1000);
    for (int i = 0; i < sampleSize; i++) {
      int value = randomData.nextBinomial(distribution.getNumberOfTrials(),
              distribution.getProbabilityOfSuccess());
      for (int j = 0; j < length; j++) {
          if (value == densityPoints[j]) {
              observedCounts[j]++;
          }
      }
    }
    TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, observedCounts, .001);
}
 

@Test
public void testNextBinomial() {
    BinomialDistributionTest testInstance = new BinomialDistributionTest();
    int[] densityPoints = testInstance.makeDensityTestPoints();
    double[] densityValues = testInstance.makeDensityTestValues();
    int sampleSize = 1000;
    int length = TestUtils.eliminateZeroMassPoints(densityPoints, densityValues);
    BinomialDistribution distribution = (BinomialDistribution) testInstance.makeDistribution();
    double[] expectedCounts = new double[length];
    long[] observedCounts = new long[length];
    for (int i = 0; i < length; i++) {
        expectedCounts[i] = sampleSize * densityValues[i];
    }
    randomData.reSeed(1000);
    for (int i = 0; i < sampleSize; i++) {
      int value = randomData.nextBinomial(distribution.getNumberOfTrials(),
              distribution.getProbabilityOfSuccess());
      for (int j = 0; j < length; j++) {
          if (value == densityPoints[j]) {
              observedCounts[j]++;
          }
      }
    }
    TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, observedCounts, .001);
}
 

@Test
public void testNextBinomial() {
    BinomialDistributionTest testInstance = new BinomialDistributionTest();
    int[] densityPoints = testInstance.makeDensityTestPoints();
    double[] densityValues = testInstance.makeDensityTestValues();
    int sampleSize = 1000;
    int length = TestUtils.eliminateZeroMassPoints(densityPoints, densityValues);
    BinomialDistribution distribution = (BinomialDistribution) testInstance.makeDistribution();
    double[] expectedCounts = new double[length];
    long[] observedCounts = new long[length];
    for (int i = 0; i < length; i++) {
        expectedCounts[i] = sampleSize * densityValues[i];
    }
    randomData.reSeed(1000);
    for (int i = 0; i < sampleSize; i++) {
      int value = randomData.nextBinomial(distribution.getNumberOfTrials(),
              distribution.getProbabilityOfSuccess());
      for (int j = 0; j < length; j++) {
          if (value == densityPoints[j]) {
              observedCounts[j]++;
          }
      }
    }
    TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, observedCounts, .001);
}
 

@Test
public void testNextBinomial() {
    BinomialDistributionTest testInstance = new BinomialDistributionTest();
    int[] densityPoints = testInstance.makeDensityTestPoints();
    double[] densityValues = testInstance.makeDensityTestValues();
    int sampleSize = 1000;
    int length = TestUtils.eliminateZeroMassPoints(densityPoints, densityValues);
    BinomialDistribution distribution = (BinomialDistribution) testInstance.makeDistribution();
    double[] expectedCounts = new double[length];
    long[] observedCounts = new long[length];
    for (int i = 0; i < length; i++) {
        expectedCounts[i] = sampleSize * densityValues[i];
    }
    randomData.reSeed(1000);
    for (int i = 0; i < sampleSize; i++) {
      int value = randomData.nextBinomial(distribution.getNumberOfTrials(),
              distribution.getProbabilityOfSuccess());
      for (int j = 0; j < length; j++) {
          if (value == densityPoints[j]) {
              observedCounts[j]++;
          }
      }
    }
    TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, observedCounts, .001);
}
 

@Test
public void testNextBinomial() {
    BinomialDistributionTest testInstance = new BinomialDistributionTest();
    int[] densityPoints = testInstance.makeDensityTestPoints();
    double[] densityValues = testInstance.makeDensityTestValues();
    int sampleSize = 1000;
    int length = TestUtils.eliminateZeroMassPoints(densityPoints, densityValues);
    BinomialDistribution distribution = (BinomialDistribution) testInstance.makeDistribution();
    double[] expectedCounts = new double[length];
    long[] observedCounts = new long[length];
    for (int i = 0; i < length; i++) {
        expectedCounts[i] = sampleSize * densityValues[i];
    }
    randomData.reSeed(1000);
    for (int i = 0; i < sampleSize; i++) {
      int value = randomData.nextBinomial(distribution.getNumberOfTrials(),
              distribution.getProbabilityOfSuccess());
      for (int j = 0; j < length; j++) {
          if (value == densityPoints[j]) {
              observedCounts[j]++;
          }
      }
    }
    TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, observedCounts, .001);
}
 

double ploidyLikelihood(double ploidy, @NotNull final WeightedPloidy weighted) {
    final String binomialKey = weighted.alleleReadCount() + ":" + weighted.totalReadCount();
    final BinomialDistribution binomialDistribution = binomialDistributionMap.computeIfAbsent(binomialKey,
            s -> new BinomialDistribution(weighted.totalReadCount(), weighted.alleleFrequency()));

    double lowerBoundAlleleReadCount = Math.max(0, ploidy - modelWidth / 2d) / weighted.ploidy() * weighted.alleleReadCount();
    int lowerBoundAlleleReadCountRounded = (int) Math.round(lowerBoundAlleleReadCount);
    double lowerBoundAddition = lowerBoundAlleleReadCountRounded + 0.5 - lowerBoundAlleleReadCount;

    double upperBoundAlleleReadCount = Math.max(0, ploidy + modelWidth / 2d) / weighted.ploidy() * weighted.alleleReadCount();
    int upperBoundAlleleReadCountRounded = (int) Math.round(upperBoundAlleleReadCount);
    double upperBoundSubtraction = upperBoundAlleleReadCountRounded + 0.5 - upperBoundAlleleReadCount;

    double rawResult =
            binomialDistribution.cumulativeProbability(upperBoundAlleleReadCountRounded) - binomialDistribution.cumulativeProbability(
                    lowerBoundAlleleReadCountRounded) + lowerBoundAddition * binomialDistribution.probability(
                    lowerBoundAlleleReadCountRounded) - upperBoundSubtraction * binomialDistribution.probability(
                    upperBoundAlleleReadCountRounded);

    return Math.round(rawResult * 100) / 100d;
}
 

@Test
public void testNextBinomial() {
    BinomialDistributionTest testInstance = new BinomialDistributionTest();
    int[] densityPoints = testInstance.makeDensityTestPoints();
    double[] densityValues = testInstance.makeDensityTestValues();
    int sampleSize = 1000;
    int length = TestUtils.eliminateZeroMassPoints(densityPoints, densityValues);
    BinomialDistribution distribution = (BinomialDistribution) testInstance.makeDistribution();
    double[] expectedCounts = new double[length];
    long[] observedCounts = new long[length];
    for (int i = 0; i < length; i++) {
        expectedCounts[i] = sampleSize * densityValues[i];
    }
    randomData.reSeed(1000);
    for (int i = 0; i < sampleSize; i++) {
      int value = randomData.nextBinomial(distribution.getNumberOfTrials(),
              distribution.getProbabilityOfSuccess());
      for (int j = 0; j < length; j++) {
          if (value == densityPoints[j]) {
              observedCounts[j]++;
          }
      }
    }
    TestUtils.assertChiSquareAccept(densityPoints, expectedCounts, observedCounts, .001);
}
 

protected static AllelicCount generateAllelicCount(final double minorFraction, final SimpleInterval position,
                                                   final RandomGenerator rng, final GammaDistribution biasGenerator, final double outlierProbability) {
    final int numReads = 100;
    final double bias = biasGenerator.sample();

    //flip a coin to decide alt minor (alt fraction = minor fraction) or ref minor (alt fraction = 1 - minor fraction)
    final double altFraction =  rng.nextDouble() < 0.5 ? minorFraction : 1 - minorFraction;

    //the probability of an alt read is the alt fraction modified by the bias or, in the case of an outlier, random
    final double pAlt = rng.nextDouble() < outlierProbability ? rng.nextDouble()
            : altFraction / (altFraction + (1 - altFraction) * bias);

    final int numAltReads = new BinomialDistribution(rng, numReads, pAlt).sample();
    final int numRefReads = numReads - numAltReads;
    return new AllelicCount(position, numAltReads, numRefReads);
}
 

/**
 * Tests slice sampling of a beta posterior = uniform prior x binomial likelihood from 1000 data points.
 * Checks that input mean and standard deviation are recovered from the posterior of the mean parameter
 * by 500 burn-in samples + 1000 samples to a relative error of 1% and 5%, respectively.
 */
@Test
public void testSliceSamplingOfBetaPosterior() {
    rng.setSeed(RANDOM_SEED);

    final int numTrials = 100;
    final double p = 0.9;
    final BinomialDistribution binomialDistribution = new BinomialDistribution(rng, numTrials, p);
    final BiFunction<Integer, Double, Double> binomialLogLikelihood =
            (d, x) -> new BinomialDistribution(null, numTrials, x).logProbability(d);
    final List<Integer> data = Ints.asList(binomialDistribution.sample(NUM_DATA_POINTS));
    final double numSuccessesTotal = data.stream().mapToDouble(x -> x).sum();
    final double alpha = 1. + numSuccessesTotal;
    final double beta = 1. + (numTrials * NUM_DATA_POINTS - numSuccessesTotal);
    final double variance = new BetaDistribution(null, alpha, beta).getNumericalVariance();

    final double xInitial = 0.5;
    final double xMin = 0.;
    final double xMax = 1.;
    final double width = 0.1;
    final int numBurnInSamples = 500;
    final int numSamples = 1500;
    final MinibatchSliceSampler<Integer> betaSampler = new MinibatchSliceSampler<>(
            rng, data, UNIFORM_LOG_PRIOR, binomialLogLikelihood,
            xMin, xMax, width, MINIBATCH_SIZE, APPROX_THRESHOLD);
    final double[] samples = Doubles.toArray(betaSampler.sample(xInitial, numSamples).subList(numBurnInSamples, numSamples));

    final double sampleMean = new Mean().evaluate(samples);
    final double sampleStandardDeviation = new StandardDeviation().evaluate(samples);
    Assert.assertEquals(relativeError(sampleMean, p), 0., 0.01);
    Assert.assertEquals(relativeError(sampleStandardDeviation, Math.sqrt(variance)), 0., 0.05);
}
 

/**
 * Calculates the theoretical sensitivity with a given Phred-scaled quality score distribution at a constant
 * depth.
 * @param depth Depth to compute sensitivity at
 * @param qualityHistogram Phred-scaled quality score histogram
 * @param logOddsThreshold Log odd threshold necessary for variant to be called
 * @param sampleSize sampleSize is the total number of simulations to run
 * @param alleleFraction the allele fraction to evaluate sensitivity at
 * @param randomSeed random number seed to use for random number generator
 * @return Theoretical sensitivity for the given arguments at a constant depth.
 */
public static double sensitivityAtConstantDepth(final int depth, final Histogram<Integer> qualityHistogram, final double logOddsThreshold, final int sampleSize, final double alleleFraction, final long randomSeed) {
    // If the depth is 0 at a particular locus, the sensitivity is trivially 0.0.
    if (depth == 0) {
        return 0.0;
    }

    final RouletteWheel qualityRW = new RouletteWheel(trimDistribution(normalizeHistogram(qualityHistogram)));
    final Random randomNumberGenerator = new Random(randomSeed);
    final RandomGenerator rg = new Well19937c(randomSeed);
    final BinomialDistribution bd = new BinomialDistribution(rg, depth, alleleFraction);

    // Calculate mean and deviation of quality score distribution to enable Gaussian sampling below
    final double averageQuality = qualityHistogram.getMean();
    final double standardDeviationQuality = qualityHistogram.getStandardDeviation();

    int calledVariants = 0;
    // Sample simulated variants, and count the number that would get called.  The ratio
    // of the number called to the total sampleSize is the sensitivity.
    for (int sample = 0; sample < sampleSize; sample++) {
        final int altDepth = bd.sample();

        final int sumOfQualities;
        if (altDepth < LARGE_NUMBER_OF_DRAWS) {
            // If the number of alt reads is "small" we draw from the actual base quality distribution.
            sumOfQualities = IntStream.range(0, altDepth).map(n -> qualityRW.draw()).sum();
        } else {
            // If the number of alt reads is "large" we draw from a Gaussian approximation of the base
            // quality distribution to speed up the code.
            sumOfQualities = drawSumOfQScores(altDepth, averageQuality, standardDeviationQuality, randomNumberGenerator.nextGaussian());
        }

        if (isCalled(depth, altDepth, (double) sumOfQualities, alleleFraction, logOddsThreshold)) {
            calledVariants++;
        }
    }
    return (double) calledVariants / sampleSize;
}
 

/**
 *
 * @param validationTumorTotalCount the total number of reads that we would like to validate.
 *                                  Must be positive or zero.
 * @param maxSignalRatioInNormal a ratio of non-ref to ref bases.  Must be within 0.0 - 1.0.
 * @return the minimum read count to be above the model noise floor.
 */
public static int calculateMinCountForSignal(int validationTumorTotalCount, double maxSignalRatioInNormal) {
    ParamUtils.isPositiveOrZero(validationTumorTotalCount, "Cannot have a negative total count.");
    ParamUtils.inRange(maxSignalRatioInNormal, 0.0, 1.0, "Cannot have have a ratio that is outside of 0.0 - 1.0.");

    final BinomialDistribution binomialDistribution = new BinomialDistribution(validationTumorTotalCount, maxSignalRatioInNormal);
    return Math.max(binomialDistribution.inverseCumulativeProbability(P_VALUE_FOR_NOISE), MINIMUM_NUM_READS_FOR_SIGNAL_COUNT);
}
 

@Override
public double calculateErrorProbability(final VariantContext vc, final Mutect2FilteringEngine filteringEngine, ReferenceContext referenceContext) {

    final int[] rpa = vc.getAttributeAsList(GATKVCFConstants.REPEATS_PER_ALLELE_KEY).stream()
            .mapToInt(o -> Integer.parseInt(String.valueOf(o))).toArray();
    if (rpa.length < 2) {
        return 0;
    }
    final String ru = vc.getAttributeAsString(GATKVCFConstants.REPEAT_UNIT_KEY, "");

    final int referenceSTRBaseCount = ru.length() * rpa[0];
    final int numPCRSlips = rpa[0] - rpa[1];
    if (referenceSTRBaseCount >= minSlippageLength && Math.abs(numPCRSlips) == 1) {
        // calculate the p-value that out of n reads we would have at least k slippage reads
        // if this p-value is small we keep the variant (reject the PCR slippage hypothesis)
        final int[] ADs = filteringEngine.sumADsOverSamples(vc, true, false);
        if (ADs == null || ADs.length < 2) {
            return 0;
        }
        final int depth = (int) MathUtils.sum(ADs);

        final int altCount = (int) MathUtils.sum(ADs) - ADs[0];
        final double logSomaticLikelihood = filteringEngine.getSomaticClusteringModel().logLikelihoodGivenSomatic(depth, altCount);

        double likelihoodGivenSlippageArtifact;
        try {
            likelihoodGivenSlippageArtifact = Beta.regularizedBeta(slippageRate, ADs[1] + 1, ADs[0] + 1);
        } catch (final MaxCountExceededException e) {
            //if the special function can't be computed, use a binomial with fixed probability
            likelihoodGivenSlippageArtifact = new BinomialDistribution(null, depth, slippageRate).probability(ADs[1]);
        }

        final double logOdds = logSomaticLikelihood - Math.log(likelihoodGivenSlippageArtifact);
        return filteringEngine.posteriorProbabilityOfError(vc, logOdds, 0);
    } else {
        return 0;
    }
}
 

@Override
public double calculateErrorProbability(final VariantContext vc, final Mutect2FilteringEngine filteringEngine, ReferenceContext referenceContext) {
    final double[] tumorLods = Mutect2FilteringEngine.getTumorLogOdds(vc);
    final int indexOfMaxTumorLod = MathUtils.maxElementIndex(tumorLods);

    final int[] tumorAlleleDepths = filteringEngine.sumADsOverSamples(vc, true, false);
    final int tumorDepth = (int) MathUtils.sum(tumorAlleleDepths);
    final int tumorAltDepth = tumorAlleleDepths[indexOfMaxTumorLod + 1];

    final int[] normalAlleleDepths = filteringEngine.sumADsOverSamples(vc, false, true);
    final int normalDepth = (int) MathUtils.sum(normalAlleleDepths);
    final int normalAltDepth = normalAlleleDepths[indexOfMaxTumorLod + 1];

    // if normal AF << tumor AF, don't filter regardless of LOD
    final double tumorAlleleFraction = (double) tumorAltDepth / tumorDepth;
    final double normalAlleleFraction = normalDepth == 0 ? 0 : (double) normalAltDepth / normalDepth;

    if (normalAlleleFraction < MIN_NORMAL_ARTIFACT_RATIO * tumorAlleleFraction)  {
        return 0.0;
    }

    final double[] normalArtifactNegativeLogOdds = MathUtils.applyToArrayInPlace(VariantContextGetters.getAttributeAsDoubleArray(vc, GATKVCFConstants.NORMAL_ARTIFACT_LOG_10_ODDS_KEY), MathUtils::log10ToLog);
    final double normalArtifactProbability = filteringEngine.posteriorProbabilityOfNormalArtifact(normalArtifactNegativeLogOdds[indexOfMaxTumorLod]);

    // the normal artifact log odds misses artifacts whose support in the normal consists entirely of low base quality reads
    // Since a lot of low-BQ reads is itself evidence of an artifact, we filter these by hand via an estimated LOD
    // that uses the average base quality of *ref* reads in the normal
    final int medianRefBaseQuality = vc.getAttributeAsIntList(GATKVCFConstants.MEDIAN_BASE_QUALITY_KEY, IMPUTED_NORMAL_BASE_QUALITY).get(0);
    final double normalPValue = 1 - new BinomialDistribution(null, normalDepth, QualityUtils.qualToErrorProb(medianRefBaseQuality))
            .cumulativeProbability(normalAltDepth - 1);

    return normalPValue < normalPileupPValueThreshold ? 1.0 : normalArtifactProbability;
}
 

/** p-value for being an artifact
 *
 * @param totalAltAlleleCount total number of alt reads
 * @param artifactAltAlleleCount alt read counts in a pair orientation that supports an artifact
 * @param biasP believed bias p value for a binomial distribution in artifact variants
 * @return p-value for this being an artifact
 */
public static double calculateArtifactPValue(final int totalAltAlleleCount, final int artifactAltAlleleCount, final double biasP) {
    ParamUtils.isPositiveOrZero(biasP, "bias parameter must be positive or zero.");
    ParamUtils.isPositiveOrZero(totalAltAlleleCount, "total alt allele count must be positive or zero.");
    ParamUtils.isPositiveOrZero(artifactAltAlleleCount, "artifact supporting alt allele count must be positive or zero.");
    ParamUtils.isPositiveOrZero(totalAltAlleleCount-artifactAltAlleleCount, "Total alt count must be same or greater than the artifact alt count.");
    return new BinomialDistribution(null, totalAltAlleleCount, biasP).cumulativeProbability(artifactAltAlleleCount);
}
 

@Test
public void testClear() throws Exception {
    final double f = 0.01;
    final int N = 100000;
    final ReadsDownsampler fd = new FractionalDownsampler(f);
    for (int i = 0; i < N; i++) {
        fd.submit(ArtificialReadUtils.createArtificialRead("100M"));
    }
    final BinomialDistribution bin = new BinomialDistribution(N, f);
    final double errorRate = 1.0 / 1000;   //we can fail this often
    Assert.assertTrue(fd.size() <= bin.inverseCumulativeProbability(1 - errorRate / 2));
    Assert.assertTrue(fd.size() >= bin.inverseCumulativeProbability(errorRate / 2));

    Assert.assertTrue(fd.hasFinalizedItems());
    Assert.assertFalse(fd.hasPendingItems());
    Assert.assertFalse(fd.requiresCoordinateSortOrder());
    Assert.assertNotNull(fd.peekFinalized());
    Assert.assertNull(fd.peekPending());
    fd.clearItems();
    Assert.assertEquals(fd.size(), 0);
    Assert.assertFalse(fd.hasFinalizedItems());
    Assert.assertFalse(fd.hasPendingItems());
    Assert.assertFalse(fd.requiresCoordinateSortOrder());
    Assert.assertNull(fd.peekFinalized());
    Assert.assertNull(fd.peekPending());

}
 

public static double expectedBAF(int averageDepth, double percent) {
    int minDepth = averageDepth / 2;
    int maxDepth = averageDepth * 3 / 2;
    double totalProbability = 0;

    for (int i = minDepth; i < maxDepth; i++) {
        final BinomialDistribution distribution = new BinomialDistribution(i, 0.5);
        double probability = 1d * distribution.inverseCumulativeProbability(percent) / i;
        totalProbability += probability;
    }

    return totalProbability / (maxDepth - minDepth);
}
 

@Override
public Object doWork(Object first, Object second) throws IOException{
  if(null == first){
    throw new IOException(String.format(Locale.ROOT,"Invalid expression %s - null found for the first value",toExpression(constructingFactory)));
  }
  if(null == second){
    throw new IOException(String.format(Locale.ROOT,"Invalid expression %s - null found for the second value",toExpression(constructingFactory)));
  }

  Number numberOfTrials = (Number)first;
  Number successProb = (Number)second;

  return new BinomialDistribution(numberOfTrials.intValue(), successProb.doubleValue());
}
 

/**
 * Test sampling from a binomial distribution.
 */
@Test
public void testBinomialSamples() {
    final int trials = 67;
    final double probabilityOfSuccess = 0.345;
    final BinomialDistribution dist = new BinomialDistribution(trials, probabilityOfSuccess);
    final double[] expected = new double[trials + 1];
    for (int i = 0; i < expected.length; i++) {
        expected[i] = dist.probability(i);
    }
    checkSamples(expected);
}
 

/**
 * Test sampling from a binomial distribution.
 */
@Test
public void testBinomialSamples() {
    final int trials = 67;
    final double probabilityOfSuccess = 0.345;
    final BinomialDistribution dist = new BinomialDistribution(null, trials, probabilityOfSuccess);
    final double[] expected = new double[trials + 1];
    for (int i = 0; i < expected.length; i++) {
        expected[i] = dist.probability(i);
    }
    checkSamples(expected, 1.0);
}
 

public static void main(String[] args) {
    BinomialDistribution bd = new BinomialDistribution(n, p);
    for (int x = 0; x <= n; x++) {
        System.out.printf("%4d%8.4f%n", x, bd.probability(x));
    }
    System.out.printf("mean = %6.4f%n", bd.getNumericalMean());
    double variance = bd.getNumericalVariance();
    double stdv = Math.sqrt(variance);
    System.out.printf("standard deviation = %6.4f%n", stdv);
}
 

static double BinomLogLik(double d, int asRef, int asAlt) {
    int total = asRef + asAlt;
    //determine likelihood here.
    BinomialDistribution binomDist = new BinomialDistribution(total, d);
    Double logDensity = binomDist.logProbability(asRef);
    return -1.0 * logDensity;
}
 

@Test(timeOut = 10_000)
public void testWeightedFairSchedulingEqualWeights()
{
    InternalResourceGroup root = new InternalResourceGroup("root", (group, export) -> {}, directExecutor()) {
        @Override
        public void triggerProcessQueuedQueries()
        {
            // No op to allow the test fine-grained control about when to trigger the next query.
        }
    };
    root.setSoftMemoryLimitBytes(DataSize.of(1, MEGABYTE).toBytes());
    root.setMaxQueuedQueries(50);
    // Start with zero capacity, so that nothing starts running until we've added all the queries
    root.setHardConcurrencyLimit(0);
    root.setSchedulingPolicy(WEIGHTED_FAIR);

    InternalResourceGroup group1 = root.getOrCreateSubGroup("1");
    group1.setSoftMemoryLimitBytes(DataSize.of(1, MEGABYTE).toBytes());
    group1.setMaxQueuedQueries(50);
    group1.setHardConcurrencyLimit(2);
    group1.setSoftConcurrencyLimit(2);
    group1.setSchedulingWeight(1);

    InternalResourceGroup group2 = root.getOrCreateSubGroup("2");
    group2.setSoftMemoryLimitBytes(DataSize.of(1, MEGABYTE).toBytes());
    group2.setMaxQueuedQueries(50);
    group2.setHardConcurrencyLimit(2);
    group2.setSoftConcurrencyLimit(2);
    group2.setSchedulingWeight(1);

    InternalResourceGroup group3 = root.getOrCreateSubGroup("3");
    group3.setSoftMemoryLimitBytes(DataSize.of(1, MEGABYTE).toBytes());
    group3.setMaxQueuedQueries(50);
    group3.setHardConcurrencyLimit(2);
    group3.setSoftConcurrencyLimit(2);
    group3.setSchedulingWeight(2);

    Set<MockManagedQueryExecution> group1Queries = fillGroupTo(group1, ImmutableSet.of(), 4);
    Set<MockManagedQueryExecution> group2Queries = fillGroupTo(group2, ImmutableSet.of(), 4);
    Set<MockManagedQueryExecution> group3Queries = fillGroupTo(group3, ImmutableSet.of(), 4);
    root.setHardConcurrencyLimit(4);

    int group1Ran = 0;
    int group2Ran = 0;
    int group3Ran = 0;
    for (int i = 0; i < 1000; i++) {
        group1Ran += completeGroupQueries(group1Queries);
        group2Ran += completeGroupQueries(group2Queries);
        group3Ran += completeGroupQueries(group3Queries);
        root.updateGroupsAndProcessQueuedQueries();
        group1Queries = fillGroupTo(group1, group1Queries, 4);
        group2Queries = fillGroupTo(group2, group2Queries, 4);
        group3Queries = fillGroupTo(group3, group3Queries, 4);
    }

    // group 3 should run approximately 2x the number of queries of 1 and 2
    BinomialDistribution binomial = new BinomialDistribution(4000, 1.0 / 4.0);
    int lowerBound = binomial.inverseCumulativeProbability(0.000001);
    int upperBound = binomial.inverseCumulativeProbability(0.999999);

    assertBetweenInclusive(group1Ran, lowerBound, upperBound);
    assertBetweenInclusive(group2Ran, lowerBound, upperBound);
    assertBetweenInclusive(group3Ran, 2 * lowerBound, 2 * upperBound);
}
 

@Override
public double calculateErrorProbability(final VariantContext vc, final Mutect2FilteringEngine filteringEngine, ReferenceContext referenceContext) {
    final double[] somaticLogOdds = Mutect2FilteringEngine.getTumorLogOdds(vc);
    final int maxLodIndex = MathUtils.maxElementIndex(somaticLogOdds);

    final Optional<double[]> normalLogOdds = vc.hasAttribute(GATKVCFConstants.NORMAL_LOG_10_ODDS_KEY) ?
            Optional.of(MathUtils.applyToArrayInPlace(VariantContextGetters.getAttributeAsDoubleArray(vc, GATKVCFConstants.NORMAL_LOG_10_ODDS_KEY), MathUtils::log10ToLog)) : Optional.empty();
    final double[] negativeLog10AlleleFrequencies = VariantContextGetters.getAttributeAsDoubleArray(vc, GATKVCFConstants.POPULATION_AF_KEY);
    final double populationAF = Math.pow(10, -negativeLog10AlleleFrequencies[maxLodIndex]);

    if (populationAF < EPSILON) {
        return 0;
    } else if (populationAF > 1 - EPSILON) {
        return 1;
    }

    // note that this includes the ref
    final int[] alleleCounts = filteringEngine.sumADsOverSamples(vc, true, false);
    final int totalCount = (int) MathUtils.sum(alleleCounts);
    if (totalCount == 0) {  // this could come up in GGA mode
        return 0;
    }
    final int refCount = alleleCounts[0];
    final int altCount = alleleCounts[maxLodIndex + 1];
    final double altAlleleFraction = filteringEngine.weightedAverageOfTumorAFs(vc)[maxLodIndex];

    final double maf = computeMinorAlleleFraction(vc, filteringEngine, alleleCounts);

    // sum of alt minor and alt major possibilities
    final double logGermlineLikelihood = NaturalLogUtils.LOG_ONE_HALF + NaturalLogUtils.logSumExp(
            new BinomialDistribution(null, totalCount, maf).logProbability(altCount),
            new BinomialDistribution(null, totalCount, 1 - maf).logProbability(altCount));

    final double logSomaticLikelihood = filteringEngine.getSomaticClusteringModel().logLikelihoodGivenSomatic(totalCount, altCount);
    // this is \chi in the docs, the correction factor for tumor likelihoods if forced to have maf or 1 - maf
    // as the allele fraction
    final double logOddsOfGermlineHetVsSomatic = logGermlineLikelihood - logSomaticLikelihood;

    // see docs -- basically the tumor likelihood for a germline hom alt is approximately equal to the somatic likelihood
    // as long as the allele fraction is high
    final double logOddsOfGermlineHomAltVsSomatic = altAlleleFraction < MIN_ALLELE_FRACTION_FOR_GERMLINE_HOM_ALT ? Double.NEGATIVE_INFINITY : 0;

    final double normalLod = normalLogOdds.isPresent() ? normalLogOdds.get()[maxLodIndex] : 0;
    // note the minus sign required because Mutect has the convention that this is log odds of allele *NOT* being in the normal
    return germlineProbability(-normalLod, logOddsOfGermlineHetVsSomatic, logOddsOfGermlineHomAltVsSomatic,
            populationAF, filteringEngine.getLogSomaticPrior(vc, maxLodIndex));
}
 

AlleleFractionSimulatedData(final SampleLocatableMetadata metadata,
                            final AlleleFractionGlobalParameters globalParameters,
                            final int numSegments,
                            final double averageHetsPerSegment,
                            final double averageDepth,
                            final RandomGenerator rng) {
    final AlleleFractionState.MinorFractions minorFractions = new AlleleFractionState.MinorFractions(numSegments);
    final List<AllelicCount> allelicCounts = new ArrayList<>();
    final List<SimpleInterval> segments = new ArrayList<>();

    final PoissonDistribution segmentLengthGenerator = makePoisson(rng, averageHetsPerSegment);
    final PoissonDistribution readDepthGenerator = makePoisson(rng, averageDepth);
    final UniformRealDistribution minorFractionGenerator = new UniformRealDistribution(rng, 0.0, 0.5);

    final double meanBias = globalParameters.getMeanBias();
    final double biasVariance = globalParameters.getBiasVariance();
    final double outlierProbability = globalParameters.getOutlierProbability();

    //translate to ApacheCommons' parametrization of the gamma distribution
    final double gammaShape = meanBias * meanBias / biasVariance;
    final double gammaScale = biasVariance / meanBias;
    final GammaDistribution biasGenerator = new GammaDistribution(rng, gammaShape, gammaScale);

    //put each segment on its own chromosome and sort in sequence-dictionary order
    final List<String> chromosomes = IntStream.range(0, numSegments)
            .mapToObj(i -> metadata.getSequenceDictionary().getSequence(i).getSequenceName())
            .collect(Collectors.toList());

    for (final String chromosome : chromosomes) {
        // calculate the range of het indices for this segment
        final int numHetsInSegment = Math.max(MIN_HETS_PER_SEGMENT, segmentLengthGenerator.sample());

        final double minorFraction = minorFractionGenerator.sample();
        minorFractions.add(minorFraction);

        //we will put all the hets in this segment/chromosome at loci 1, 2, 3 etc
        segments.add(new SimpleInterval(chromosome, 1, numHetsInSegment));
        for (int het = 1; het < numHetsInSegment + 1; het++) {
            final double bias = biasGenerator.sample();

            //flip a coin to decide alt minor (alt fraction = minor fraction) or ref minor (alt fraction = 1 - minor fraction)
            final boolean isAltMinor = rng.nextDouble() < 0.5;
            final double altFraction =  isAltMinor ? minorFraction : 1 - minorFraction;

            //the probability of an alt read is the alt fraction modified by the bias or, in the case of an outlier, random
            final double pAlt;
            if (rng.nextDouble() < outlierProbability) {
                pAlt = rng.nextDouble();
            } else {
                pAlt = altFraction / (altFraction + (1 - altFraction) * bias);
            }

            final int numReads = readDepthGenerator.sample();
            final int numAltReads = new BinomialDistribution(rng, numReads, pAlt).sample();
            final int numRefReads = numReads - numAltReads;
            allelicCounts.add(new AllelicCount(new SimpleInterval(chromosome, het, het), numRefReads, numAltReads));
        }
    }

    data = new AlleleFractionSegmentedData(
            new AllelicCountCollection(metadata, allelicCounts),
            new SimpleIntervalCollection(metadata, segments));
    trueState = new AlleleFractionState(meanBias, biasVariance, outlierProbability, minorFractions);
}
 

@Test
public void testMLUpdate() throws Exception {
  Path tempDir = getTempDir();
  Path dataDir = tempDir.resolve("data");
  Map<String,Object> overlayConfig = new HashMap<>();
  overlayConfig.put("oryx.batch.update-class", MockMLUpdate.class.getName());
  ConfigUtils.set(overlayConfig, "oryx.batch.storage.data-dir", dataDir);
  ConfigUtils.set(overlayConfig, "oryx.batch.storage.model-dir", tempDir.resolve("model"));
  overlayConfig.put("oryx.batch.streaming.generation-interval-sec", GEN_INTERVAL_SEC);
  overlayConfig.put("oryx.ml.eval.test-fraction", TEST_FRACTION);
  overlayConfig.put("oryx.ml.eval.threshold", DATA_TO_WRITE / 2); // Should easily pass threshold
  Config config = ConfigUtils.overlayOn(overlayConfig, getConfig());

  startMessaging();

  List<Integer> trainCounts = MockMLUpdate.getResetTrainCounts();
  List<Integer> testCounts = MockMLUpdate.getResetTestCounts();

  startServerProduceConsumeTopics(config, DATA_TO_WRITE, WRITE_INTERVAL_MSEC);

  // If lists are unequal at this point, there must have been an empty test set
  // which yielded no call to evaluate(). Fill in the blank
  while (trainCounts.size() > testCounts.size()) {
    testCounts.add(0);
  }

  log.info("trainCounts = {}", trainCounts);
  log.info("testCounts = {}", testCounts);

  checkOutputData(dataDir, DATA_TO_WRITE);
  checkIntervals(trainCounts.size(), DATA_TO_WRITE, WRITE_INTERVAL_MSEC, GEN_INTERVAL_SEC);

  assertEquals(testCounts.size(), trainCounts.size());

  RandomGenerator random = RandomManager.getRandom();
  int lastTotalTrainCount = 0;
  int lastTestCount = 0;
  for (int i = 0; i < testCounts.size(); i++) {
    int totalTrainCount = trainCounts.get(i);
    int testCount = testCounts.get(i);
    int newTrainInGen = totalTrainCount - (lastTotalTrainCount + lastTestCount);
    if (newTrainInGen == 0) {
      continue;
    }
    lastTotalTrainCount = totalTrainCount;
    lastTestCount = testCount;
    int totalNew = testCount + newTrainInGen;

    IntegerDistribution dist = new BinomialDistribution(random, totalNew, TEST_FRACTION);
    checkDiscreteProbability(testCount, dist);
  }

}
 

@Test
public void testBernoulliLogProb() {
    Nd4j.getRandom().setSeed(12345);

    int inputSize = 4;
    int[] mbs = new int[] {1, 2, 5};

    Random r = new Random(12345);

    for (boolean average : new boolean[] {true, false}) {
        for (int minibatch : mbs) {

            INDArray x = Nd4j.zeros(minibatch, inputSize);
            for (int i = 0; i < minibatch; i++) {
                for (int j = 0; j < inputSize; j++) {
                    x.putScalar(i, j, r.nextInt(2));
                }
            }

            INDArray distributionParams = Nd4j.rand(minibatch, inputSize).muli(2).subi(1); //i.e., pre-sigmoid prob
            INDArray prob = Transforms.sigmoid(distributionParams, true);

            ReconstructionDistribution dist = new BernoulliReconstructionDistribution(Activation.SIGMOID);

            double negLogProb = dist.negLogProbability(x, distributionParams, average);

            INDArray exampleNegLogProb = dist.exampleNegLogProbability(x, distributionParams);
            assertArrayEquals(new long[] {minibatch, 1}, exampleNegLogProb.shape());

            //Calculate the same thing, but using Apache Commons math

            double logProbSum = 0.0;
            for (int i = 0; i < minibatch; i++) {
                double exampleSum = 0.0;
                for (int j = 0; j < inputSize; j++) {
                    double p = prob.getDouble(i, j);

                    BinomialDistribution binomial = new BinomialDistribution(1, p); //Bernoulli is a special case of binomial

                    double xVal = x.getDouble(i, j);
                    double thisLogProb = binomial.logProbability((int) xVal);
                    logProbSum += thisLogProb;
                    exampleSum += thisLogProb;
                }
                assertEquals(-exampleNegLogProb.getDouble(i), exampleSum, 1e-6);
            }

            double expNegLogProb;
            if (average) {
                expNegLogProb = -logProbSum / minibatch;
            } else {
                expNegLogProb = -logProbSum;
            }

            //                System.out.println(x);

            //                System.out.println(expNegLogProb + "\t" + logProb + "\t" + (logProb / expNegLogProb));
            assertEquals(expNegLogProb, negLogProb, 1e-6);

            //Also: check random sampling...
            int count = minibatch * inputSize;
            INDArray arr = Nd4j.linspace(-3, 3, count, Nd4j.dataType()).reshape(minibatch, inputSize);
            INDArray sampleMean = dist.generateAtMean(arr);
            INDArray sampleRandom = dist.generateRandom(arr);

            for (int i = 0; i < minibatch; i++) {
                for (int j = 0; j < inputSize; j++) {
                    double d1 = sampleMean.getDouble(i, j);
                    double d2 = sampleRandom.getDouble(i, j);
                    assertTrue(d1 >= 0.0 || d1 <= 1.0); //Mean value - probability... could do 0 or 1 (based on most likely) but that isn't very useful...
                    assertTrue(d2 == 0.0 || d2 == 1.0);
                }
            }
        }
    }
}
 

/**
 * Returns the <i>observed significance level</i>, or
 * <a href="http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">p-value</a>,
 * associated with a <a href="http://en.wikipedia.org/wiki/Binomial_test"> Binomial test</a>.
 * <p>
 * The number returned is the smallest significance level at which one can reject the null hypothesis.
 * The form of the hypothesis depends on {@code alternativeHypothesis}.</p>
 * <p>
 * The p-Value represents the likelihood of getting a result at least as extreme as the sample,
 * given the provided {@code probability} of success on a single trial. For single-sided tests,
 * this value can be directly derived from the Binomial distribution. For the two-sided test,
 * the implementation works as follows: we start by looking at the most extreme cases
 * (0 success and n success where n is the number of trials from the sample) and determine their likelihood.
 * The lower value is added to the p-Value (if both values are equal, both are added). Then we continue with
 * the next extreme value, until we added the value for the actual observed sample.</p>
 * <p>
 * <strong>Preconditions</strong>:
 * <ul>
 * <li>Number of trials must be &ge; 0.</li>
 * <li>Number of successes must be &ge; 0.</li>
 * <li>Number of successes must be &le; number of trials.</li>
 * <li>Probability must be &ge; 0 and &le; 1.</li>
 * </ul></p>
 *
 * @param numberOfTrials number of trials performed
 * @param numberOfSuccesses number of successes observed
 * @param probability assumed probability of a single trial under the null hypothesis
 * @param alternativeHypothesis type of hypothesis being evaluated (one- or two-sided)
 * @return p-value
 * @throws NotPositiveException if {@code numberOfTrials} or {@code numberOfSuccesses} is negative
 * @throws OutOfRangeException if {@code probability} is not between 0 and 1
 * @throws MathIllegalArgumentException if {@code numberOfTrials} &lt; {@code numberOfSuccesses} or
 * if {@code alternateHypothesis} is null.
 * @see AlternativeHypothesis
 */
public double binomialTest(int numberOfTrials, int numberOfSuccesses, double probability,
                           AlternativeHypothesis alternativeHypothesis) {
    if (numberOfTrials < 0) {
        throw new NotPositiveException(numberOfTrials);
    }
    if (numberOfSuccesses < 0) {
        throw new NotPositiveException(numberOfSuccesses);
    }
    if (probability < 0 || probability > 1) {
        throw new OutOfRangeException(probability, 0, 1);
    }
    if (numberOfTrials < numberOfSuccesses) {
        throw new MathIllegalArgumentException(
            LocalizedFormats.BINOMIAL_INVALID_PARAMETERS_ORDER,
            numberOfTrials, numberOfSuccesses);
    }
    if (alternativeHypothesis == null) {
        throw new NullArgumentException();
    }

    // pass a null rng to avoid unneeded overhead as we will not sample from this distribution
    final BinomialDistribution distribution = new BinomialDistribution(null, numberOfTrials, probability);
    switch (alternativeHypothesis) {
    case GREATER_THAN:
        return 1 - distribution.cumulativeProbability(numberOfSuccesses - 1);
    case LESS_THAN:
        return distribution.cumulativeProbability(numberOfSuccesses);
    case TWO_SIDED:
        int criticalValueLow = 0;
        int criticalValueHigh = numberOfTrials;
        double pTotal = 0;

        while (true) {
            double pLow = distribution.probability(criticalValueLow);
            double pHigh = distribution.probability(criticalValueHigh);

            if (pLow == pHigh) {
                pTotal += 2 * pLow;
                criticalValueLow++;
                criticalValueHigh--;
            } else if (pLow < pHigh) {
                pTotal += pLow;
                criticalValueLow++;
            } else {
                pTotal += pHigh;
                criticalValueHigh--;
            }

            if (criticalValueLow > numberOfSuccesses || criticalValueHigh < numberOfSuccesses) {
                break;
            }
        }
        return pTotal;
    default:
        throw new MathInternalError(LocalizedFormats. OUT_OF_RANGE_SIMPLE, alternativeHypothesis,
                  AlternativeHypothesis.TWO_SIDED, AlternativeHypothesis.LESS_THAN);
    }
}
 

/**
 * Returns the <i>observed significance level</i>, or
 * <a
 * href="http://www.cas.lancs.ac.uk/glossary_v1.1/hyptest.html#pvalue">p-value</a>,
 * associated with a <a href="http://en.wikipedia.org/wiki/Binomial_test">
 * Binomial test</a>.
 * <p>
 * The number returned is the smallest significance level at which one can
 * reject the null hypothesis. The form of the hypothesis depends on
 * {@code alternativeHypothesis}.</p>
 * <p>
 * The p-Value represents the likelihood of getting a result at least as
 * extreme as the sample, given the provided {@code probability} of success
 * on a single trial. For single-sided tests, this value can be directly
 * derived from the Binomial distribution. For the two-sided test, the
 * implementation works as follows: we start by looking at the most extreme
 * cases (0 success and n success where n is the number of trials from the
 * sample) and determine their likelihood. The lower value is added to the
 * p-Value (if both values are equal, both are added). Then we continue with
 * the next extreme value, until we added the value for the actual observed
 * sample.</p>
 * <p>
 * <strong>Preconditions</strong>:
 * <ul>
 * <li>Number of trials must be &ge; 0.</li>
 * <li>Number of successes must be &ge; 0.</li>
 * <li>Number of successes must be &le; number of trials.</li>
 * <li>Probability must be &ge; 0 and &le; 1.</li>
 * </ul></p>
 *
 * @param numberOfTrials number of trials performed
 * @param numberOfSuccesses number of successes observed
 * @param probability assumed probability of a single trial under the null
 * hypothesis
 * @param alternativeHypothesis type of hypothesis being evaluated (one- or
 * two-sided)
 * @return p-value
 * @throws NotPositiveException if {@code numberOfTrials} or
 * {@code numberOfSuccesses} is negative
 * @throws OutOfRangeException if {@code probability} is not between 0 and 1
 * @throws MathIllegalArgumentException if {@code numberOfTrials} &lt;
 * {@code numberOfSuccesses} or if {@code alternateHypothesis} is null.
 * @see AlternativeHypothesis
 */
public double binomialTest(int numberOfTrials, int numberOfSuccesses, double probability,
		AlternativeHypothesis alternativeHypothesis) {
	if (numberOfTrials < 0) {
		throw new NotPositiveException(numberOfTrials);
	}
	if (numberOfSuccesses < 0) {
		throw new NotPositiveException(numberOfSuccesses);
	}
	if (probability < 0 || probability > 1) {
		throw new OutOfRangeException(probability, 0, 1);
	}
	if (numberOfTrials < numberOfSuccesses) {
		throw new MathIllegalArgumentException(
				LocalizedFormats.BINOMIAL_INVALID_PARAMETERS_ORDER,
				numberOfTrials, numberOfSuccesses);
	}
	if (alternativeHypothesis == null) {
		throw new NullArgumentException();
	}

	final BinomialDistribution distribution = new BinomialDistribution(numberOfTrials, probability);
	switch (alternativeHypothesis) {
		case GREATER_THAN:
			return 1 - distribution.cumulativeProbability(numberOfSuccesses - 1);
		case LESS_THAN:
			return distribution.cumulativeProbability(numberOfSuccesses);
		case TWO_SIDED:
			int criticalValueLow = 0;
			int criticalValueHigh = numberOfTrials;
			double pTotal = 0;

			while (true) {
				double pLow = distribution.probability(criticalValueLow);
				double pHigh = distribution.probability(criticalValueHigh);

				if (pLow == pHigh) {
					pTotal += 2 * pLow;
					criticalValueLow++;
					criticalValueHigh--;
				} else if (pLow < pHigh) {
					pTotal += pLow;
					criticalValueLow++;
				} else {
					pTotal += pHigh;
					criticalValueHigh--;
				}

				if (criticalValueLow > numberOfSuccesses || criticalValueHigh < numberOfSuccesses) {
					break;
				}
			}
			return pTotal;
		default:
			throw new MathInternalError(LocalizedFormats.OUT_OF_RANGE_SIMPLE, alternativeHypothesis,
					AlternativeHypothesis.TWO_SIDED, AlternativeHypothesis.LESS_THAN);
	}
}