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

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

@operator (
		value = "gamma_density",
		can_be_const = false,
		category = { IOperatorCategory.RANDOM },
		concept = { IConcept.RANDOM })
@doc (
		value = "gamma_density(x,shape,scale) returns the probability density function (PDF) at the specified point x "
				+ "of the Gamma distribution with the given shape and scale.",
		examples = { @example ( value = "gamma_density(1,9,0.5)", equals = "0.731", test = false) },
		see = { "binomial", "gauss_rnd", "lognormal_rnd", "poisson", "rnd", "skew_gauss", "truncated_gauss",
				"weibull_rnd", "weibull_density", "lognormal_density" })
@no_test (Reason.IMPOSSIBLE_TO_TEST)
public static Double OpGammaDist(final IScope scope, final Double x, 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.density(x);
}
 

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

@Override
public void prepare(Configurator configuration, DataModule model, DataSpace space) {
    super.prepare(configuration, model, space);
    userMean = configuration.getFloat("recommender.recommender.user.mu", 0F);
    userBeta = configuration.getFloat("recommender.recommender.user.beta", 1F);
    userWishart = configuration.getFloat("recommender.recommender.user.wishart.scale", 1F);

    itemMean = configuration.getFloat("recommender.recommender.item.mu", 0F);
    itemBeta = configuration.getFloat("recommender.recommender.item.beta", 1F);
    itemWishart = configuration.getFloat("recommender.recommender.item.wishart.scale", 1F);

    rateSigma = configuration.getFloat("recommender.recommender.rating.sigma", 2F);

    gibbsIterations = configuration.getInteger("recommender.recommender.iterations.gibbs", 1);

    userMatrixes = new DenseMatrix[epocheSize - 1];
    itemMatrixes = new DenseMatrix[epocheSize - 1];

    normalDistribution = new QuantityProbability(JDKRandomGenerator.class, factorSize, NormalDistribution.class, 0D, 1D);
    userGammaDistributions = new QuantityProbability[factorSize];
    itemGammaDistributions = new QuantityProbability[factorSize];
    for (int index = 0; index < factorSize; index++) {
        userGammaDistributions[index] = new QuantityProbability(JDKRandomGenerator.class, index, GammaDistribution.class, (userSize + factorSize - (index + 1D)) / 2D, 2D);
        itemGammaDistributions[index] = new QuantityProbability(JDKRandomGenerator.class, index, GammaDistribution.class, (itemSize + factorSize - (index + 1D)) / 2D, 2D);
    }
}
 

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

@Override
public Double map(Row row) throws Exception {
    double label = (Double) row.getField(numFeature);
    double weight = (Double) row.getField(numFeature + 1);
    double pred = (Double) row.getField(numFeature + 3);

    GammaDistribution distribution = new GammaDistribution(1.0 / disp, 1 / (pred * disp));
    return weight * Math.log(distribution.density(label));
}
 

@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 shape = (Number)first;
  Number scale = (Number)second;

  return new GammaDistribution(shape.doubleValue(), scale.doubleValue());
}
 

@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 OnlineLDAModel(int K, float alpha, float eta, long D, double tau0, double kappa,
        double delta) {
    super(K);

    if (tau0 < 0.d) {
        throw new IllegalArgumentException("tau0 MUST be positive: " + tau0);
    }
    if (kappa <= 0.5 || 1.d < kappa) {
        throw new IllegalArgumentException("kappa MUST be in (0.5, 1.0]: " + kappa);
    }

    this._alpha = alpha;
    this._eta = eta;
    this._D = D;
    this._tau0 = tau0;
    this._kappa = kappa;
    this._delta = delta;

    this._isAutoD = (_D <= 0L);

    // initialize a random number generator
    this._gd = new GammaDistribution(SHAPE, SCALE);
    _gd.reseedRandomGenerator(1001);

    // initialize the parameters
    this._lambda = new HashMap<String, float[]>(100);
}
 

protected static AllelicCountCollection generateCounts(final List<Double> minorAlleleFractionSequence,
                                                       final List<SimpleInterval> positions,
                                                       final RandomGenerator rng,
                                                       final AlleleFractionGlobalParameters trueParams) {
    //translate to ApacheCommons' parametrization of the gamma distribution
    final GammaDistribution biasGenerator = getGammaDistribution(trueParams, rng);
    final double outlierProbability = trueParams.getOutlierProbability();

    final AllelicCountCollection counts = new AllelicCountCollection();

    for (int n = 0; n < minorAlleleFractionSequence.size(); n++) {
        counts.add(generateAllelicCount(minorAlleleFractionSequence.get(n), positions.get(n), rng, biasGenerator, outlierProbability));
    }
    return counts;
}
 

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

@Test
public void testBasics() throws IOException {
    // this sampler has four variables
    // g1 is gamma distributed with alpha = 0.2, beta = 0.2
    // v1 is unit normal
    // v2 is normal with mean = 0, sd = 2
    // v3 is gamma-normal with dof=2, mean = 0.
    SchemaSampler s = new SchemaSampler(Resources.asCharSource(Resources.getResource("schema015.json"), Charsets.UTF_8).read());

    TDigest tdG1 = new AVLTreeDigest(500);
    TDigest tdG2 = new AVLTreeDigest(500);
    TDigest td1 = new AVLTreeDigest(500);
    TDigest td2 = new AVLTreeDigest(500);
    TDigest td3 = new AVLTreeDigest(500);

    double x1 = 0;
    double x2 = 0;
    double x3 = 0;

    for (int i = 0; i < 1000000; i++) {
        JsonNode r = s.sample();
        tdG1.add(r.get("g1").asDouble());
        tdG2.add(r.get("g2").asDouble());

        double step1 = r.get("v1").get("step").asDouble();
        td1.add(step1);
        x1 += step1;
        assertEquals(x1, r.get("v1").get("value").asDouble(), 0);
        assertEquals(x1, r.get("v1-bare").asDouble(), 0);

        double step2 = r.get("v2").get("step").asDouble();
        td2.add(step2);
        x2 += step2;
        assertEquals(x2, r.get("v2").get("value").asDouble(), 0);

        double step3 = r.get("v3").get("step").asDouble();
        td3.add(step3);
        x3 += step3;
        assertEquals(x3, r.get("v3").get("value").asDouble(), 0);
    }

    // now compare against reference distributions to test accuracy of the observed step distributions
    NormalDistribution normalDistribution = new NormalDistribution();
    GammaDistribution gd1 = new GammaDistribution(0.2, 5);
    GammaDistribution gd2 = new GammaDistribution(1, 1);
    TDistribution tDistribution = new TDistribution(2);
    for (double q : new double[]{0.001, 0.01, 0.1, 0.2, 0.5, 0.8, 0.9, 0.99, 0.99}) {
        double uG1 = gd1.cumulativeProbability(tdG1.quantile(q));
        assertEquals(q, uG1, (1 - q) * q * 10e-2);

        double uG2 = gd2.cumulativeProbability(tdG2.quantile(q));
        assertEquals(q, uG2, (1 - q) * q * 10e-2);

        double u1 = normalDistribution.cumulativeProbability(td1.quantile(q));
        assertEquals(q, u1, (1 - q) * q * 10e-2);

        double u2 = normalDistribution.cumulativeProbability(td2.quantile(q) / 2);
        assertEquals(q, u2, (1 - q) * q * 10e-2);

        double u3 = tDistribution.cumulativeProbability(td3.quantile(q));
        assertEquals(q, u3, (1 - q) * q * 10e-2);
    }
}
 

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

protected static GammaDistribution getGammaDistribution(final AlleleFractionGlobalParameters trueParams, final RandomGenerator rng) {
    final double gammaShape = trueParams.getMeanBias() * trueParams.getMeanBias() / trueParams.getBiasVariance();
    final double gammaScale = trueParams.getBiasVariance() / trueParams.getMeanBias();
    return new GammaDistribution(rng, gammaShape, gammaScale);
}
 

public QAlphaUnivariateFunction(final double[] effectivePhis) {
    numEffectivePhis = effectivePhis.length;
    sumLogEffectivePhis = Arrays.stream(effectivePhis).map(Math::log).sum();
    gammaDistribution = new GammaDistribution(null, GAMMA_SHAPE, GAMMA_SCALE);
}
 

public Pair<Double, Double> sampleWithNormalGamma(final N state) {
	double tau = new GammaDistribution(this.alpha.get(state), this.beta.get(state)).sample();
	double muNew = new NormalDistribution(this.mu.get(state), 1 / (this.lambda.get(state) * tau)).sample();
	return new Pair<>(muNew, tau);
}
 

private static Tuple4<DenseMatrix, DenseMatrix, Long, Long> onlineCorpusUpdate(
    List<Vector> data, DenseMatrix lambda, DenseMatrix alpha, DenseMatrix gammad,
    int vocabularySize, int numTopic, double subSamplingRate) {

    boolean isRandGamma = gammad == null;
    //wordTopicStat is the word topic probability information.
    DenseMatrix wordTopicStat = DenseMatrix.zeros(numTopic, vocabularySize);
    DenseMatrix logPhatPart = new DenseMatrix(numTopic, 1);
    DenseMatrix expELogBeta = LdaUtil.expDirichletExpectation(lambda).transpose();
    long nonEmptyWordCount = 0;
    long nonEmptyDocCount = 0;
    //the online corpus update stage can update the model in two way.
    //if the document order is determined, then it will update in the order.
    //or it will choose documents randomly.
    RandomDataGenerator random = new RandomDataGenerator();
    GammaDistribution distribution = new GammaDistribution(100, 100);
    int dataSize = data.size();
    boolean sampled = false;
    for (int j = 0; j < dataSize; ++j) {
        //if the subSamplingRate is too small and no doc is sampled in one iteration, then will randomly
        //choose one doc to update the model.
        double rate = random.nextUniform(0, 1);
        int index = -1;
        if (rate < subSamplingRate) {
            index = j;
            sampled = true;
        }
        if (j + 1 == dataSize && !sampled) {
            index = random.nextInt(0, dataSize - 1);
        }
        if (index != -1) {
            Vector vec = data.get(index);
            SparseVector sv = (SparseVector) vec;
            sv.setSize(vocabularySize);
            sv.removeZeroValues();
            for (int i = 0; i < sv.numberOfValues(); i++) {
                nonEmptyWordCount += sv.getValues()[i];
            }
            if (isRandGamma) {
                if (gammad == null) {
                    gammad = new DenseMatrix(numTopic, 1);
                }
                for (int i = 0; i < numTopic; i++) {
                    gammad.set(i, 0, distribution.inverseCumulativeProbability(random.nextUniform(0, 1)));
                }
            }
            Tuple2<DenseMatrix, DenseMatrix> topicDistributionTuple =
                LdaUtil.getTopicDistributionMethod(sv, expELogBeta, alpha, gammad, numTopic);

            for (int i = 0; i < sv.getIndices().length; i++) {
                for (int k = 0; k < numTopic; k++) {
                    wordTopicStat.add(k, sv.getIndices()[i], topicDistributionTuple.f1.get(k, i));
                }
            }
            gammad = topicDistributionTuple.f0;
            DenseMatrix deGammad = LdaUtil.dirichletExpectationVec(topicDistributionTuple.f0);
            for (int k = 0; k < numTopic; k++) {
                logPhatPart.add(k, 0, deGammad.get(k, 0));
            }
            nonEmptyDocCount++;
        }
    }
    return new Tuple4<>(wordTopicStat, logPhatPart, nonEmptyWordCount, nonEmptyDocCount);
}
 

/**
 * <p>Generates a random value from the
 * {@link org.apache.commons.math3.distribution.GammaDistribution Gamma Distribution}.</p>
 *
 * <p>This implementation uses the following algorithms: </p>
 *
 * <p>For 0 < shape < 1: <br/>
 * Ahrens, J. H. and Dieter, U., <i>Computer methods for
 * sampling from gamma, beta, Poisson and binomial distributions.</i>
 * Computing, 12, 223-246, 1974.</p>
 *
 * <p>For shape >= 1: <br/>
 * Marsaglia and Tsang, <i>A Simple Method for Generating
 * Gamma Variables.</i> ACM Transactions on Mathematical Software,
 * Volume 26 Issue 3, September, 2000.</p>
 *
 * @param shape the median of the Gamma distribution
 * @param scale the scale parameter of the Gamma distribution
 * @return random value sampled from the Gamma(shape, scale) distribution
 * @throws NotStrictlyPositiveException if {@code shape <= 0} or
 * {@code scale <= 0}.
 */
public double nextGamma(double shape, double scale) throws NotStrictlyPositiveException {
    return new GammaDistribution(getRandomGenerator(),shape, scale,
            GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY).sample();
}
 

/**
 * <p>Generates a random value from the
 * {@link org.apache.commons.math3.distribution.GammaDistribution Gamma Distribution}.</p>
 *
 * <p>This implementation uses the following algorithms: </p>
 *
 * <p>For 0 < shape < 1: <br/>
 * Ahrens, J. H. and Dieter, U., <i>Computer methods for
 * sampling from gamma, beta, Poisson and binomial distributions.</i>
 * Computing, 12, 223-246, 1974.</p>
 *
 * <p>For shape >= 1: <br/>
 * Marsaglia and Tsang, <i>A Simple Method for Generating
 * Gamma Variables.</i> ACM Transactions on Mathematical Software,
 * Volume 26 Issue 3, September, 2000.</p>
 *
 * @param shape the median of the Gamma distribution
 * @param scale the scale parameter of the Gamma distribution
 * @return random value sampled from the Gamma(shape, scale) distribution
 * @throws NotStrictlyPositiveException if {@code shape <= 0} or
 * {@code scale <= 0}.
 */
public double nextGamma(double shape, double scale) throws NotStrictlyPositiveException {
    return new GammaDistribution(getRan(),shape, scale,
            GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY).sample();
}
 

/**
 * <p>Generates a random value from the
 * {@link org.apache.commons.math3.distribution.GammaDistribution Gamma Distribution}.</p>
 *
 * <p>This implementation uses the following algorithms: </p>
 *
 * <p>For 0 < shape < 1: <br/>
 * Ahrens, J. H. and Dieter, U., <i>Computer methods for
 * sampling from gamma, beta, Poisson and binomial distributions.</i>
 * Computing, 12, 223-246, 1974.</p>
 *
 * <p>For shape >= 1: <br/>
 * Marsaglia and Tsang, <i>A Simple Method for Generating
 * Gamma Variables.</i> ACM Transactions on Mathematical Software,
 * Volume 26 Issue 3, September, 2000.</p>
 *
 * @param shape the median of the Gamma distribution
 * @param scale the scale parameter of the Gamma distribution
 * @return random value sampled from the Gamma(shape, scale) distribution
 * @throws NotStrictlyPositiveException if {@code shape <= 0} or
 * {@code scale <= 0}.
 */
public double nextGamma(double shape, double scale) throws NotStrictlyPositiveException {
    return new GammaDistribution(getRandomGenerator(),shape, scale,
            GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY).sample();
}
 

/**
 * <p>Generates a random value from the
 * {@link org.apache.commons.math3.distribution.GammaDistribution Gamma Distribution}.</p>
 *
 * <p>This implementation uses the following algorithms: </p>
 *
 * <p>For 0 < shape < 1: <br/>
 * Ahrens, J. H. and Dieter, U., <i>Computer methods for
 * sampling from gamma, beta, Poisson and binomial distributions.</i>
 * Computing, 12, 223-246, 1974.</p>
 *
 * <p>For shape >= 1: <br/>
 * Marsaglia and Tsang, <i>A Simple Method for Generating
 * Gamma Variables.</i> ACM Transactions on Mathematical Software,
 * Volume 26 Issue 3, September, 2000.</p>
 *
 * @param shape the median of the Gamma distribution
 * @param scale the scale parameter of the Gamma distribution
 * @return random value sampled from the Gamma(shape, scale) distribution
 * @throws NotStrictlyPositiveException if {@code shape <= 0} or
 * {@code scale <= 0}.
 */
public double nextGamma(double shape, double scale) throws NotStrictlyPositiveException {
    return new GammaDistribution(getRandomGenerator(),shape, scale,
            GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY).sample();
}
 

/**
 * <p>Generates a random value from the
 * {@link org.apache.commons.math3.distribution.GammaDistribution Gamma Distribution}.</p>
 *
 * <p>This implementation uses the following algorithms: </p>
 *
 * <p>For 0 < shape < 1: <br/>
 * Ahrens, J. H. and Dieter, U., <i>Computer methods for
 * sampling from gamma, beta, Poisson and binomial distributions.</i>
 * Computing, 12, 223-246, 1974.</p>
 *
 * <p>For shape >= 1: <br/>
 * Marsaglia and Tsang, <i>A Simple Method for Generating
 * Gamma Variables.</i> ACM Transactions on Mathematical Software,
 * Volume 26 Issue 3, September, 2000.</p>
 *
 * @param shape the median of the Gamma distribution
 * @param scale the scale parameter of the Gamma distribution
 * @return random value sampled from the Gamma(shape, scale) distribution
 * @throws NotStrictlyPositiveException if {@code shape <= 0} or
 * {@code scale <= 0}.
 */
public double nextGamma(double shape, double scale) throws NotStrictlyPositiveException {
    return new GammaDistribution(getRandomGenerator(),shape, scale,
            GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY).sample();
}
 

/**
 * <p>Generates a random value from the
 * {@link org.apache.commons.math3.distribution.GammaDistribution Gamma Distribution}.</p>
 *
 * <p>This implementation uses the following algorithms: </p>
 *
 * <p>For 0 < shape < 1: <br/>
 * Ahrens, J. H. and Dieter, U., <i>Computer methods for
 * sampling from gamma, beta, Poisson and binomial distributions.</i>
 * Computing, 12, 223-246, 1974.</p>
 *
 * <p>For shape >= 1: <br/>
 * Marsaglia and Tsang, <i>A Simple Method for Generating
 * Gamma Variables.</i> ACM Transactions on Mathematical Software,
 * Volume 26 Issue 3, September, 2000.</p>
 *
 * @param shape the median of the Gamma distribution
 * @param scale the scale parameter of the Gamma distribution
 * @return random value sampled from the Gamma(shape, scale) distribution
 * @throws NotStrictlyPositiveException if {@code shape <= 0} or
 * {@code scale <= 0}.
 */
public double nextGamma(double shape, double scale) throws NotStrictlyPositiveException {
    return new GammaDistribution(getRandomGenerator(),shape, scale,
            GammaDistribution.DEFAULT_INVERSE_ABSOLUTE_ACCURACY).sample();
}