I was working on something like this but I'll stop since you've got it. I think this is the way to go. If you are going to be making several calls to it you can save the threshold list by something like
poissonthreshlist =: $0 poissonran=: 4 : 0"0 if. (#poissonthreshlist) ~: n =. 5*x do. poissonthreshlist =: +/\(^-x)**/\1,x%}.i.>:n end. poissonthreshlist I. y [EMAIL PROTECTED] 0 ) As a final tweak you could look to see if there are any values of (#threshlist) returned - there usually won't be - and replace them with values drawn from the tail of the distribution. Henry Rich > -----Original Message----- > From: [EMAIL PROTECTED] > [mailto:[EMAIL PROTECTED] On Behalf Of Robert Cyr > Sent: Friday, August 08, 2008 4:22 PM > To: Programming forum > Subject: Re: [Jprogramming] general Gamma distribution > > 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 ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
