Java Code Examples for org.apache.commons.math3.distribution.GammaDistribution#sample()

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);
}
 

@operator (
		value = "gamma_rnd",
		can_be_const = false,
		category = { IOperatorCategory.RANDOM },
		concept = { IConcept.RANDOM })
@doc (
		value = "returns a random value from a gamma distribution with specified values of the shape and scale parameters",
		examples = { @example (
				value = "gamma_rnd(9,0.5)",
				equals = "0.731",
				test = false) },
		see = { "binomial", "gauss_rnd", "lognormal_rnd", "poisson", "rnd", "skew_gauss", "truncated_gauss",
				"weibull_rnd", "gamma_trunc_rnd" })
@no_test (Reason.IMPOSSIBLE_TO_TEST)
public static Double OpGammaDist(final IScope scope, final Double shape, final Double scale)
		throws GamaRuntimeException {
	final GammaDistribution dist = new GammaDistribution(scope.getRandom().getGenerator(), shape, scale,
			GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
	return dist.sample();
}
 

public static double[] staticDraw(double[] alpha) {
    double dist[] = new double[alpha.length];
    
    // For each dimension, draw a sample from Gamma(mp_i, 1).
    for (int i = 0; i < dist.length; i++) {
        GammaDistribution gammaDist = new GammaDistribution(rng, alpha[i], 1, 
                GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY);
        dist[i] = gammaDist.sample();
        if (dist[i] <= 0) {
            dist[i] = EPSILON;
        }
    }
    
    // Normalize the distribution.
    Multinomials.normalizeProps(dist);

    return dist;
}
 

@Nonnull
public static float[] newRandomFloatArray(@Nonnegative final int size,
        @Nonnull final GammaDistribution gd) {
    final float[] ret = new float[size];
    for (int i = 0; i < size; i++) {
        ret[i] = (float) gd.sample();
    }
    return ret;
}
 

public AlleleFractionSimulatedData(final double averageHetsPerSegment, final int numSegments,
        final double averageDepth, final double biasMean, final double biasVariance, final double outlierProbability) {
    rng.setSeed(RANDOM_SEED);
    this.numSegments = numSegments;
    final AlleleFractionState.MinorFractions minorFractions = new AlleleFractionState.MinorFractions(numSegments);
    final List<AllelicCount> alleleCounts = 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);

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

    //put each segment on its own chromosome and sort by lexicographical order
    final List<String> chromosomes = IntStream.range(0, numSegments).mapToObj(Integer::toString).collect(Collectors.toList());
    Collections.sort(chromosomes);

    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 + 1));
        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) {
                truePhases.add(AlleleFractionIndicator.OUTLIER);
                pAlt = rng.nextDouble();
            } else {
                truePhases.add(isAltMinor ? AlleleFractionIndicator.ALT_MINOR : AlleleFractionIndicator.REF_MINOR);
                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;
            alleleCounts.add(new AllelicCount(new SimpleInterval(chromosome, het, het), numRefReads, numAltReads));
        }
    }

    final Genome genome = new Genome(TRIVIAL_TARGETS, alleleCounts);
    segmentedGenome = new SegmentedGenome(segments, genome);
    trueState = new AlleleFractionState(biasMean, biasVariance, outlierProbability, minorFractions);
}
 

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);
}