(Random rand)
| 98 | |
| 99 | |
| 100 | private double sampleOne(Random rand) |
| 101 | { |
| 102 | //From http://www.johndcook.com/blog/2010/06/14/generating-poisson-random-values/ |
| 103 | double c = 0.767 - 3.36/lambda; |
| 104 | double beta = PI/sqrt(3.0*lambda); |
| 105 | double alpha = beta*lambda; |
| 106 | double k = log(c) - lambda - log(beta); |
| 107 | |
| 108 | while(true) |
| 109 | { |
| 110 | double u = rand.nextDouble(); |
| 111 | double x = (alpha - log((1.0 - u) / u)) / beta; |
| 112 | double n = floor(x + 0.5); |
| 113 | if (n < 0) |
| 114 | continue; |
| 115 | double v = rand.nextDouble(); |
| 116 | double y = alpha - beta * x; |
| 117 | // double lhs = y + log(v/(1.0 + exp(y))^2); |
| 118 | //simplify right part as log(v)-2 log(e^y+1) |
| 119 | // double lhs = y + log(v/pow(1.0 + exp(y), 2)); |
| 120 | double lhs = y + log(v) - 2 * log(exp(y) + 1); |
| 121 | // double rhs = k + n*log(lambda) - log(n!); |
| 122 | double rhs = k + n * log(lambda) - SpecialMath.lnGamma(n + 1); |
| 123 | if (lhs <= rhs) |
| 124 | return n; |
| 125 | } |
| 126 | } |
| 127 | |
| 128 | @Override |
| 129 | public double[] sample(int numSamples, Random rand) |
no test coverage detected