java cmos_Java PoissonDistribution.sample方法代码示例

import org.apache.commons.math3.distribution.PoissonDistribution; //导入方法依赖的package包/类

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 allelicCounts = new ArrayList<>();

final List 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 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);

}