/** * {@inheritDoc} * <p> * <strong>Algorithm Description</strong>: * <ul><li> For small means, uses simulation of a Poisson process * using Uniform deviates, as described * <a href="http://irmi.epfl.ch/cmos/Pmmi/interactive/rng7.htm"> here.</a> * The Poisson process (and hence value returned) is bounded by 1000 * mean.</li> * * <li> For large means, uses the rejection algorithm described in <br/> * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i> * <strong>Computing</strong> vol. 26 pp. 197-207.</li></ul></p> * @throws NotStrictlyPositiveException if {@code len <= 0} */ public long nextPoisson(double mean) throws NotStrictlyPositiveException { return new PoissonDistribution(getRandomGenerator(), mean, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS).sample(); }
double sample(double u) { if (u < limit) { List<WeightedThing<Integer>> steps = Lists.newArrayList(); limit = 1; int i = 0; while (u / 20 < limit) { double pdf = pd.probability(i); limit -= pdf; steps.add(new WeightedThing<>(i, pdf)); i++; } steps.add(new WeightedThing<>(steps.size(), limit)); partial = new Multinomial<>(steps); } return partial.sample(u); } }
/** * {@inheritDoc} * * For mean parameter {@code p}, the variance is {@code p}. */ public double getNumericalVariance() { return getMean(); }
/** * Create a poisson sampler which can sample elements with replacement. * * @param fraction The expected count of each element. * @param seed Random number generator seed for internal PoissonDistribution. */ public PoissonSampler(double fraction, long seed) { Preconditions.checkArgument(fraction >= 0, "fraction should be positive."); this.fraction = fraction; if (this.fraction > 0) { this.poissonDistribution = new PoissonDistribution(fraction); this.poissonDistribution.reseedRandomGenerator(seed); } this.random = new XORShiftRandom(seed); }
/** * Create a poisson sampler which can sample elements with replacement. * * @param fraction The expected count of each element. */ public PoissonSampler(double fraction) { Preconditions.checkArgument(fraction >= 0, "fraction should be non-negative."); this.fraction = fraction; if (this.fraction > 0) { this.poissonDistribution = new PoissonDistribution(fraction); } this.random = new XORShiftRandom(); }
static public double probability(PoissonDistribution poissonDistribution, Number x){ org.apache.commons.math3.distribution.PoissonDistribution distribution = new org.apache.commons.math3.distribution.PoissonDistribution(null, poissonDistribution.getMean(), org.apache.commons.math3.distribution.PoissonDistribution.DEFAULT_EPSILON, org.apache.commons.math3.distribution.PoissonDistribution.DEFAULT_MAX_ITERATIONS); x = (Number)TypeUtil.cast(DataType.INTEGER, x); return distribution.probability(x.intValue()); }
@Test public void logProbMatchesKnownLogProb() { DoubleVertex mu = ConstantVertex.of(new double[]{5, 7}); PoissonVertex poissonVertex = new PoissonVertex(mu); double actualLogProb = poissonVertex.logPmf(new int[]{39, 49}); PoissonDistribution distribution1 = new PoissonDistribution(5); PoissonDistribution distribution2 = new PoissonDistribution(7); double expectedLogProb = distribution1.logProbability(39) + distribution2.logProbability(49); assertEquals(expectedLogProb, actualLogProb, 1e-5); } }
/** * @param param * mean * @return Poisson distribution */ protected PoissonDistribution getPoissonDistribution(double param) { if (poisson == null || poisson.getMean() != param) { poisson = new PoissonDistribution(param); } return poisson; }
private void samplingProcess() { if (fraction <= THRESHOLD) { double u = Math.max(random.nextDouble(), EPSILON); int gap = (int) (Math.log(u) / -fraction); skipGapElements(gap); if (input.hasNext()) { currentElement = input.next(); currentCount = poisson_ge1(fraction); } } else { while (input.hasNext()) { currentElement = input.next(); currentCount = poissonDistribution.sample(); if (currentCount > 0) { break; } } } } };
/** {@inheritDoc} */ public double probability(int x) { final double logProbability = logProbability(x); return logProbability == Double.NEGATIVE_INFINITY ? 0 : FastMath.exp(logProbability); }
@Override public final void compute() { if (input[0].isDefined() && input[1].isDefined() && input[2].isDefined()) { double param = a.getDouble(); int val = (int) Math.round(b.getDouble()); try { PoissonDistribution dist = getPoissonDistribution(param); if (isCumulative.getBoolean()) { num.setValue(dist.cumulativeProbability(val)); // P(X <= // val) } else { num.setValue(dist.probability(val)); // P(X = val) } } catch (Exception e) { Log.debug(e.getMessage()); num.setUndefined(); } } else { num.setUndefined(); } }
public PoissonSampler(double lambda) { limit = 1; gen = RandomUtils.getRandom(); pd = new PoissonDistribution(gen.getRandomGenerator(), lambda, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS); }
/** * Create a poisson sampler which can sample elements with replacement. * * @param fraction The expected count of each element. * @param seed Random number generator seed for internal PoissonDistribution. */ public Poisson(double fraction, long seed) { Preconditions.checkArgument(fraction >= 0, "fraction should be positive."); this.fraction = fraction; if (this.fraction > 0) { this.poissonDistribution = new PoissonDistribution(fraction); this.poissonDistribution.reseedRandomGenerator(seed); } this.random = new XORShiftRNG(seed); }
/** * Checks whether the observations conform to a Poisson process with the * specified intensity. Uses a chi square test with the specified confidence. * The null hypothesis is that the observations are the result of a poisson * process. * @param observations * @param intensity * @param confidence * @return <code>true</code> if the observations */ static boolean isPoissonProcess(Frequency observations, double intensity, double length, double confidence) { final PoissonDistribution pd = new PoissonDistribution(length * intensity); final Iterator<?> it = observations.valuesIterator(); final long[] observed = new long[observations.getUniqueCount()]; final double[] expected = new double[observations.getUniqueCount()]; int index = 0; while (it.hasNext()) { final Long l = (Long) it.next(); observed[index] = observations.getCount(l); expected[index] = pd.probability(l.intValue()) * observations.getSumFreq(); if (expected[index] == 0) { return false; } index++; } final double chi = TestUtils.chiSquareTest(expected, observed); return !(chi < confidence); }
@Test public void logProbForValuesGreaterThanTwenty() { double mu = 25.0; PoissonVertex poissonVertex = new PoissonVertex(mu); double logProb19 = poissonVertex.logPmf(19); double logProb20 = poissonVertex.logPmf(20); double logProb100 = poissonVertex.logPmf(100); PoissonDistribution distribution = new PoissonDistribution(25); assertThat(logProb19, closeTo(distribution.logProbability(19), 1e-6)); assertThat(logProb20, closeTo(distribution.logProbability(20), 1e-6)); assertThat(logProb100, closeTo(distribution.logProbability(100), 1e-6)); }
private void samplingProcess() { if (fraction <= THRESHOLD) { double u = Math.max(random.nextDouble(), EPSILON); int gap = (int) (Math.log(u) / -fraction); skipGapElements(gap); if (input.hasNext()) { currentElement = input.next(); currentCount = poisson_ge1(fraction); } } else { while (input.hasNext()) { currentElement = input.next(); currentCount = poissonDistribution.sample(); if (currentCount > 0) { break; } } } } };
/** {@inheritDoc} */ public double probability(int x) { final double logProbability = logProbability(x); return logProbability == Double.NEGATIVE_INFINITY ? 0 : Math.exp(logProbability); }
/** {@inheritDoc} */ @Override public Stream<UpstreamEntry> transform(Stream<UpstreamEntry> upstream) { PoissonDistribution poisson = new PoissonDistribution( new Well19937c(seed), subsampleRatio, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS); return upstream.sequential().flatMap(en -> Stream.generate(() -> en).limit(poisson.sample())); } }
private static void checkDistribution(Sampler<Double> pd, double alpha) { int[] count = new int[(int) Math.max(10, 5 * alpha)]; for (int i = 0; i < 10000; i++) { count[pd.sample().intValue()]++; } IntegerDistribution ref = new PoissonDistribution(RandomUtils.getRandom().getRandomGenerator(), alpha, PoissonDistribution.DEFAULT_EPSILON, PoissonDistribution.DEFAULT_MAX_ITERATIONS); for (int i = 0; i < count.length; i++) { assertEquals(ref.probability(i), count[i] / 10000.0, 2.0e-2); } } }
/** * Create a poisson sampler which can sample elements with replacement. * * @param fraction The expected count of each element. * @param seed Random number generator seed for internal PoissonDistribution. */ public PoissonSampler(double fraction, long seed) { Preconditions.checkArgument(fraction >= 0, "fraction should be positive."); this.fraction = fraction; if (this.fraction > 0) { this.poissonDistribution = new PoissonDistribution(fraction); this.poissonDistribution.reseedRandomGenerator(seed); } this.random = new XORShiftRandom(seed); }