Robert,
I believe the recursive relation to which you refer
is this.
prob[X=x+1] = lambda*(prob[X=x])%(x+1)
In Fishman's book the first algorithm is like yours,
but does not seem to guess at 'n=. 5*x'
0 Set omega to exp[-lambda] and make A and B double
precision
1 IF lambda<0, omega<:0, or omega>:1, go to 7
2 X=. 0
A=. omega
B=. A
3 Generate U from ?0
4 if. (U<:A) +. (A>0.999999) +. (B<0.000001)
do.
X
return.
end.
5 X=. >:X
6 B=. B*lambda%X
A=. A+B
go to 4
7 Print error message stop
On Fri, 8 Aug 2008, Robert Cyr wrote:
+ I have googled a few methods for generating Poisson random variable. They
+ typically base their calculation of the distribution function of a Poisson
+ using a recursive process , cumulate the result until a randomly chosen
+ number is found. They then return the position.
+
+ This is compatible with the poissondist technique which one can easily
+ recognize as the standard Poisson distribution function. Applying the above
+ technique leads to:
+
+ poissonran=: 4 : 0
+ n=. 5*x
+ (+/\(^-x)**/\1,x%}.i.>:n)I.?y$0
+ )
+
+ It's really fast and accurate for integer values as long as the limit is
+ chosen sufficiently large to ensure that the last of the cumulated values is
+ 1.
+
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
