@@ -58,8 +58,9 @@ impl Sample<u64> for Poisson {
58
58
impl IndependentSample < u64 > for Poisson {
59
59
fn ind_sample < R : Rng > ( & self , rng : & mut R ) -> u64 {
60
60
// using the algorithm from Numerical Recipes in C
61
+
62
+ // for low expected values use the Knuth method
61
63
if self . lambda < 12.0 {
62
- // for lambda below the threshold, use the Knuth method
63
64
let mut result = 0 ;
64
65
let mut p = 1.0 ;
65
66
while p > self . exp_lambda {
@@ -68,23 +69,38 @@ impl IndependentSample<u64> for Poisson {
68
69
}
69
70
result - 1
70
71
}
72
+ // high expected values - rejection method
71
73
else {
72
- // rejection method
74
+ // some cached values
73
75
let tmp = ( 2.0 * self . lambda ) . sqrt ( ) ;
74
76
let mut int_result: u64 ;
77
+
75
78
loop {
76
79
let mut result;
77
80
let mut comp_dev;
78
- // first, generate an easier deviate greater than 0
81
+
82
+ // we use the lorentzian distribution as the comparison distribution
83
+ // f(x) ~ 1/(1+x/^2)
79
84
loop {
80
- comp_dev = ( PI * rng. gen :: < f64 > ( ) ) . tan ( ) ; // comparison deviate
85
+ // draw from the lorentzian distribution
86
+ comp_dev = ( PI * rng. gen :: < f64 > ( ) ) . tan ( ) ;
87
+ // shift the peak of the comparison ditribution
81
88
result = tmp* comp_dev + self . lambda ;
89
+ // repeat the drawing until we are in the range of possible values
82
90
if result >= 0.0 { break ; }
83
91
}
84
- // now result is a random variable greater than 0 with Lorentzian distribution
92
+ // now the result is a random variable greater than 0 with Lorentzian distribution
93
+ // the result should be an integer value
85
94
result = result. floor ( ) ;
86
95
int_result = result as u64 ;
96
+
97
+ // this is the ratio of the Poisson distribution to the comparison distribution
98
+ // the magic value scales the distribution function to a range of approximately 0-1
99
+ // since it is not exact, we multiply the ratio by 0.9 to avoid ratios greater than 1
100
+ // this doesn't change the resulting distribution, only increases the rate of failed drawings
87
101
let check = 0.9 * ( 1.0 + comp_dev* comp_dev) * ( result* self . log_lambda - log_gamma ( 1.0 + result) - self . magic_val ) . exp ( ) ;
102
+
103
+ // check with uniform random value - if below the threshold, we are within the target distribution
88
104
if rng. gen :: < f64 > ( ) <= check { break ; }
89
105
}
90
106
int_result
@@ -110,7 +126,7 @@ mod test {
110
126
println ! ( "Poisson average: {}" , avg) ;
111
127
assert ! ( ( avg - 10.0 ) . abs( ) < 0.5 ) ; // not 100% certain, but probable enough
112
128
}
113
-
129
+
114
130
#[ test]
115
131
#[ should_panic]
116
132
#[ cfg_attr( target_env = "msvc" , ignore) ]
0 commit comments