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. I am still studying Raul Miller's proposal. Robert Cyr On Fri, Aug 8, 2008 at 3:37 PM, Oleg Kobchenko <[EMAIL PROTECTED]> wrote: > > From: Oleg Kobchenko <[EMAIL PROTECTED]> > > > > > From: Brian Schott > > > > > > 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. > > > > dstat poissonrand 100 10000 > > sample size: 10000 > > minimum: 64 > > maximum: 141 > > median: 100 > > mean: 100.006 > > std devn: 10.0345 > > skewness: 0.12287 > > kurtosis: 3.03206 > > dstat <.100+10*normalrand 10000 > > sample size: 10000 > > minimum: 61 > > maximum: 137 > > median: 99 > > mean: 99.3547 > > std devn: 10.0235 > > skewness: 0.0298283 > > kurtosis: 3.0363 > > > > 3 ts '$<.100+10*normalrand 10000' > > 0.00326381 1.31347e6 > > 3 ts '$poissonrand 100 10000' > > 0.513158 690880 > > You could compare the two samples, visually it > is strikingly close. > > load 'stats plot' > > plots=: 4 : 0 > pd 'reset' > pd (~. ; #/.~) /:~ poissonrand x , y > pd (~. ; #/.~) /:~ <.0.5+x+(%:x)*normalrand y > pd 'show' > ) > > 100 plots 10000 > > BTW, is there a library for J to do parametric or non-parametric > sample comparisons? > > The closest I got is: > > N regression & (/:~) P > Var. Coeff. S.E. t > 0 1.18791 0.05245 22.65 > 1 0.99323 0.00052 1892.61 > > Source D.F. S.S. M.S. F > Regression 1 978915.01515 978915.01515 3581955.71 > Error 9998 2732.35995 0.27329 > Total 9999 981647.37510 > > S.E. of estimate 0.52277 > Corr. coeff. squared 0.99722 > > > > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
