[Jprogramming] general Gamma distribution

[Henry Rich]

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 ?@$ 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: programming-bounces at jsoftware.com 
> [mailto:programming-bounces at jsoftware.com] 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 
> <olegykj at yahoo.com> wrote:
> 
> > > From: Oleg Kobchenko <olegykj at yahoo.com>
> > >
> > > > 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