001/* 002 * Licensed to the Apache Software Foundation (ASF) under one or more 003 * contributor license agreements. See the NOTICE file distributed with 004 * this work for additional information regarding copyright ownership. 005 * The ASF licenses this file to You under the Apache License, Version 2.0 006 * (the "License"); you may not use this file except in compliance with 007 * the License. You may obtain a copy of the License at 008 * 009 * http://www.apache.org/licenses/LICENSE-2.0 010 * 011 * Unless required by applicable law or agreed to in writing, software 012 * distributed under the License is distributed on an "AS IS" BASIS, 013 * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. 014 * See the License for the specific language governing permissions and 015 * limitations under the License. 016 */ 017package org.apache.commons.rng.sampling.distribution; 018 019import org.apache.commons.rng.UniformRandomProvider; 020 021/** 022 * Sampler for the <a href="http://mathworld.wolfram.com/PoissonDistribution.html">Poisson distribution</a>. 023 * 024 * <ul> 025 * <li> 026 * For small means, a Poisson process is simulated using uniform deviates, as 027 * described <a href="http://mathaa.epfl.ch/cours/PMMI2001/interactive/rng7.htm">here</a>. 028 * The Poisson process (and hence, the returned value) is bounded by 1000 * mean. 029 * </li> 030 * <li> 031 * For large means, we use the rejection algorithm described in 032 * <blockquote> 033 * Devroye, Luc. (1981).<i>The Computer Generation of Poisson Random Variables</i><br> 034 * <strong>Computing</strong> vol. 26 pp. 197-207. 035 * </blockquote> 036 * </li> 037 * </ul> 038 */ 039public class PoissonSampler 040 extends SamplerBase 041 implements DiscreteSampler { 042 /** Value for switching sampling algorithm. */ 043 private static final double PIVOT = 40; 044 /** Mean of the distribution. */ 045 private final double mean; 046 /** Exponential. */ 047 private final ContinuousSampler exponential; 048 /** Gaussian. */ 049 private final ContinuousSampler gaussian; 050 /** {@code log(n!)}. */ 051 private final InternalUtils.FactorialLog factorialLog; 052 053 /** 054 * @param rng Generator of uniformly distributed random numbers. 055 * @param mean Mean. 056 * @throws IllegalArgumentException if {@code mean <= 0}. 057 */ 058 public PoissonSampler(UniformRandomProvider rng, 059 double mean) { 060 super(rng); 061 if (mean <= 0) { 062 throw new IllegalArgumentException(mean + " <= " + 0); 063 } 064 065 this.mean = mean; 066 067 gaussian = new BoxMullerGaussianSampler(rng, 0, 1); 068 exponential = new AhrensDieterExponentialSampler(rng, 1); 069 factorialLog = mean < PIVOT ? 070 null : // Not used. 071 InternalUtils.FactorialLog.create().withCache((int) Math.min(mean, 2 * PIVOT)); 072 } 073 074 /** {@inheritDoc} */ 075 @Override 076 public int sample() { 077 return (int) Math.min(nextPoisson(mean), Integer.MAX_VALUE); 078 } 079 080 /** {@inheritDoc} */ 081 @Override 082 public String toString() { 083 return "Poisson deviate [" + super.toString() + "]"; 084 } 085 086 /** 087 * @param meanPoisson Mean. 088 * @return the next sample. 089 */ 090 private long nextPoisson(double meanPoisson) { 091 if (meanPoisson < PIVOT) { 092 double p = Math.exp(-meanPoisson); 093 long n = 0; 094 double r = 1; 095 096 while (n < 1000 * meanPoisson) { 097 r *= nextDouble(); 098 if (r >= p) { 099 n++; 100 } else { 101 break; 102 } 103 } 104 return n; 105 } else { 106 final double lambda = Math.floor(meanPoisson); 107 final double lambdaFractional = meanPoisson - lambda; 108 final double logLambda = Math.log(lambda); 109 final double logLambdaFactorial = factorialLog((int) lambda); 110 final long y2 = lambdaFractional < Double.MIN_VALUE ? 0 : nextPoisson(lambdaFractional); 111 final double delta = Math.sqrt(lambda * Math.log(32 * lambda / Math.PI + 1)); 112 final double halfDelta = delta / 2; 113 final double twolpd = 2 * lambda + delta; 114 final double a1 = Math.sqrt(Math.PI * twolpd) * Math.exp(1 / (8 * lambda)); 115 final double a2 = (twolpd / delta) * Math.exp(-delta * (1 + delta) / twolpd); 116 final double aSum = a1 + a2 + 1; 117 final double p1 = a1 / aSum; 118 final double p2 = a2 / aSum; 119 final double c1 = 1 / (8 * lambda); 120 121 double x = 0; 122 double y = 0; 123 double v = 0; 124 int a = 0; 125 double t = 0; 126 double qr = 0; 127 double qa = 0; 128 while (true) { 129 final double u = nextDouble(); 130 if (u <= p1) { 131 final double n = gaussian.sample(); 132 x = n * Math.sqrt(lambda + halfDelta) - 0.5d; 133 if (x > delta || x < -lambda) { 134 continue; 135 } 136 y = x < 0 ? Math.floor(x) : Math.ceil(x); 137 final double e = exponential.sample(); 138 v = -e - 0.5 * n * n + c1; 139 } else { 140 if (u > p1 + p2) { 141 y = lambda; 142 break; 143 } else { 144 x = delta + (twolpd / delta) * exponential.sample(); 145 y = Math.ceil(x); 146 v = -exponential.sample() - delta * (x + 1) / twolpd; 147 } 148 } 149 a = x < 0 ? 1 : 0; 150 t = y * (y + 1) / (2 * lambda); 151 if (v < -t && a == 0) { 152 y = lambda + y; 153 break; 154 } 155 qr = t * ((2 * y + 1) / (6 * lambda) - 1); 156 qa = qr - (t * t) / (3 * (lambda + a * (y + 1))); 157 if (v < qa) { 158 y = lambda + y; 159 break; 160 } 161 if (v > qr) { 162 continue; 163 } 164 if (v < y * logLambda - factorialLog((int) (y + lambda)) + logLambdaFactorial) { 165 y = lambda + y; 166 break; 167 } 168 } 169 return y2 + (long) y; 170 } 171 } 172 173 /** 174 * Compute the natural logarithm of the factorial of {@code n}. 175 * 176 * @param n Argument. 177 * @return {@code log(n!)} 178 * @throws IllegalArgumentException if {@code n < 0}. 179 */ 180 private double factorialLog(int n) { 181 return factorialLog.value(n); 182 } 183}