Example usage for org.apache.commons.math3.distribution BinomialDistribution sample

Introduction

In this page you can find the example usage for org.apache.commons.math3.distribution BinomialDistribution sample.

Prototype

public int sample() 

Source Link

Document

The default implementation uses the <a href="http://en.wikipedia.org/wiki/Inverse_transform_sampling"> inversion method</a>.

Usage

From source file:math.test.TestDimpleRandom.java

/**
 * Measures speed of Apache vs Colt generators
 *//*from w  ww .  java 2 s  .  c  om*/
public static void main(String[] args) {
    RandomGenerator apacheGenerator = new MersenneTwister(42);
    cern.jet.random.engine.RandomEngine cernGenerator = new cern.jet.random.engine.MersenneTwister();

    final int N = 100000;
    long start, end;

    start = System.nanoTime();
    for (int i = N; --i >= 0;) {
        apacheGenerator.nextDouble();
    }
    end = System.nanoTime();
    long apacheTime = end - start;

    start = System.nanoTime();
    for (int i = N; --i >= 0;) {
        cernGenerator.nextDouble();
    }
    end = System.nanoTime();

    long cernTime = end - start;

    System.out.format("MersenneTwister.nextDouble() x %d apache/cern %f\n", N, (double) apacheTime / cernTime);

    start = System.nanoTime();
    for (int i = N; --i >= 0;) {
        BetaDistribution apacheBeta = new BetaDistribution(apacheGenerator, 1.0, 1.0);
        apacheBeta.sample();
    }
    end = System.nanoTime();
    apacheTime = end - start;

    cern.jet.random.Beta cernBeta = new cern.jet.random.Beta(1, 1, cernGenerator);

    start = System.nanoTime();
    for (int i = N; --i >= 0;) {
        cernBeta.nextDouble();
    }
    end = System.nanoTime();
    apacheTime = end - start;

    System.out.format("Beta x %d apache/cern %f\n", N, (double) apacheTime / cernTime);

    start = System.nanoTime();
    for (int i = N; --i >= 0;) {
        GammaDistribution apacheGamma = new GammaDistribution(apacheGenerator, 1.0, 1.0);
        apacheGamma.sample();
    }
    end = System.nanoTime();
    apacheTime = end - start;

    cern.jet.random.Gamma cernGamma = new cern.jet.random.Gamma(1, 1, cernGenerator);

    start = System.nanoTime();
    for (int i = N; --i >= 0;) {
        cernGamma.nextDouble();
    }
    end = System.nanoTime();
    apacheTime = end - start;

    System.out.format("Gamma x %d apache/cern %f\n", N, (double) apacheTime / cernTime);

    start = System.nanoTime();
    for (int i = N; --i >= 0;) {
        BinomialDistribution apacheBinomial = new BinomialDistribution(apacheGenerator, 1, .5);
        apacheBinomial.sample();
    }
    end = System.nanoTime();
    apacheTime = end - start;

    cern.jet.random.Binomial cernBinomial = new cern.jet.random.Binomial(1, .5, cernGenerator);

    start = System.nanoTime();
    for (int i = N; --i >= 0;) {
        cernBinomial.nextInt();
    }
    end = System.nanoTime();
    apacheTime = end - start;

    System.out.format("Binomial x %d apache/cern %f\n", N, (double) apacheTime / cernTime);
}

---

From source file:edu.oregonstate.eecs.mcplan.domains.inventory.InventorySimulator.java

public static void applyDynamics(final InventoryState s) {
    s.r = 0;/*from  ww w  .  jav  a2s  . c o m*/
    for (int i = 0; i < s.problem.Nproducts; ++i) {
        // Pending orders arrive
        final BinomialDistribution f = new BinomialDistribution(s.orders[i], s.problem.delivery_probability);
        final int arrivals = f.sample();
        //         System.out.println( "\tarrivals[" + i + "] = " + arrivals );
        s.orders[i] -= arrivals;
        s.inventory[i] += arrivals;

        // Demand is satisfied
        final int supplied = Math.min(s.demand[i], s.inventory[i]);
        s.inventory[i] -= supplied;
        s.r += s.problem.price[i] * supplied;

        // Excess inventory is wasted
        s.inventory[i] = Math.min(s.inventory[i], s.problem.max_inventory);

        // Warehouse cost
        s.r -= s.inventory[i] * s.problem.warehouse_cost;

        // Demand changes
        //         s.demand[i] = s.rng.nextInt( s.problem.max_demand + 1 );
    }

    // Demand changes
    final int[] new_demand = s.problem.sampleNextDemand(s);
    Fn.memcpy(s.demand, new_demand);
}

---

From source file:com.analog.lyric.math.DimpleRandom.java

/**
 * Returns sample from beta distribution with specified alpha and beta parameters.
 * @since 0.08//from   ww w.j a va 2 s.  c  o m
 */
public int nextBinomial(int n, double p) {
    // randBinomial doesn't accept zero N value or 1 or 0 p value
    if (n <= 0)
        return 0;
    else if (p <= 0)
        return 0;
    else if (p >= 1)
        return n;

    BinomialDistribution randBinomial = _randBinomial;

    if (randBinomial.getNumberOfTrials() != n || randBinomial.getProbabilityOfSuccess() != p) {
        _randBinomial = randBinomial = new BinomialDistribution(_randGenerator, n, p);
    }

    return randBinomial.sample();
}

---

From source file:org.asoem.greyfish.utils.collect.BitString.java

/**
 * Create a random bit string of given {@code length} where each bit is set with probability {@code p}, and not set
 * with probability {@code 1-p}.//from   w ww . j av a2s  .  com
 *
 * @param length the length of the bit string
 * @param rng    the random number generator to use
 * @param p      the probability for each bit in the new bit string to hold the value 1
 * @return a new bit string
 */
public static BitString random(final int length, final RandomGenerator rng, final double p) {
    checkNotNull(rng);
    checkArgument(p >= 0 && p <= 1);
    checkArgument(length >= 0);

    if (length == 0) {
        return emptyBitSequence();
    }

    if (p == 0.5) {
        return random(length, rng); // faster
    }

    final int n;
    if (p == 0) {
        n = 0;
    } else if (p == 1) {
        n = length;
    } else {
        final BinomialDistribution binomialDistribution = new BinomialDistribution(rng, length, p);
        n = binomialDistribution.sample();
    }
    assert n >= 0 && n <= length : n;

    if (n == 0) {
        return zeros(length);
    } else if (n == length) {
        return ones(length);
    }

    final ContiguousSet<Integer> indexRange = ContiguousSet.create(Range.closedOpen(0, length),
            DiscreteDomain.integers());
    final Iterable<Integer> uniqueIndexSample = Samplings.random(rng).withoutReplacement().sample(indexRange,
            n);

    if ((double) n / length < 1.0 / 32) { // < 1 bit per word?
        return new IndexSetString(ImmutableSortedSet.copyOf(uniqueIndexSample), length);
    } else {
        final BitSet bs = new BitSet(length);
        for (Integer index : uniqueIndexSample) {
            bs.set(index, true);
        }
        return new BitSetString(bs, length);
    }
}

---

From source file:org.asoem.greyfish.utils.collect.BitStringTest.java

public static Map<Integer, Integer> expectedFrequencies(final int length, final double p, final int sampleSize,
        final RandomGenerator rng) {
    final Map<Integer, Integer> frequencies = Maps.newHashMap();
    final BinomialDistribution binomialDistribution = new BinomialDistribution(rng, length, p);
    for (int i = 0; i < sampleSize; i++) {
        final int sample = binomialDistribution.sample();
        Integer integer = frequencies.get(sample);
        if (integer == null) {
            frequencies.put(sample, 0);// w ww .  j a  va  2 s  .com
            integer = 0;
        }
        frequencies.put(sample, integer + 1);
    }
    return frequencies;
}

---

From source file:picard.analysis.TheoreticalSensitivity.java

/**
 * Calculates the theoretical sensitivity with a given Phred-scaled quality score distribution at a constant
 * depth.// w w  w .ja  va 2  s .  co  m
 * @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) {
    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;
}

---