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

The following examples show how to use org.apache.commons.math3.distribution.GammaDistribution. 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: Dirichlet.java From pacaya with Apache License 2.0

6 votes

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

Example #2

Source File: Random.java From gama with GNU General Public License v3.0

6 votes

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

Example #3

Source File: AlleleFractionSegmenterUnitTest.java From gatk-protected with BSD 3-Clause "New" or "Revised" License

6 votes

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

Example #4

Source File: BPMFModel.java From jstarcraft-rns with Apache License 2.0

6 votes

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

Example #5

Source File: Random.java From gama with GNU General Public License v3.0

6 votes

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

Example #6

Source File: GlmUtil.java From Alink with Apache License 2.0

5 votes

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

Example #7

Source File: GammaDistributionEvaluator.java From lucene-solr with Apache License 2.0

5 votes

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

Example #8

Source File: ArrayUtils.java From incubator-hivemall with Apache License 2.0

5 votes

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

Example #9

Source File: OnlineLDAModel.java From incubator-hivemall with Apache License 2.0

5 votes

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

Example #10

Source File: AlleleFractionSegmenterUnitTest.java From gatk-protected with BSD 3-Clause "New" or "Revised" License

5 votes

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

Example #11

Source File: AlleleFractionSimulatedData.java From gatk-protected with BSD 3-Clause "New" or "Revised" License

4 votes

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

Example #12

Source File: RandomWalkSamplerTest.java From log-synth with Apache License 2.0

4 votes

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

Example #13

Source File: AlleleFractionSimulatedData.java From gatk with BSD 3-Clause "New" or "Revised" License

4 votes

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

Example #14

Source File: AlleleFractionSegmenterUnitTest.java From gatk-protected with BSD 3-Clause "New" or "Revised" License

4 votes

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

Example #15

Source File: CNLOHCaller.java From gatk-protected with BSD 3-Clause "New" or "Revised" License

4 votes

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

Example #16

Source File: DNGPolicy.java From AILibs with GNU Affero General Public License v3.0

4 votes

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

Example #17

Source File: OnlineCorpusStep.java From Alink with Apache License 2.0

4 votes

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

Example #18