Robert,
The poisson distribution is almost identical to the
normal distribution when its parameter is 15 or larger, so
you may want to use the normal distribution. More
specifically, according to George Fishman, as the mean,
lambda, increases the distribution of
(X-lambda)/sqrt(lambda) converges to the standard normal
distribution (ie mean 0, std dev 1); if Y is from standard
normal then
X=max(0,integer [lambda+Y*sqrt(lambda)+0.5])
A key feature would be to generate the poisson
cumulative probabilities only one time as you did in your
revised approach, whether you use the Poisson or the normal
approximation. If the mean of the Poisson does not vary,
then why not generate all 1 000 000 at once?
It seems to me that you could generate 1000000
Poisson and uniform variates all at one time (in an array of
shape 500 2000?) and then focus on keeping track of the
(ending) inventory level (i) in your main process using J's
strength of array processing. This is the challenging part,
imo, especially avoiding the for. for iteration.
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm