[Jprogramming] general Gamma distribution

[Roger Hui]

tpdf2=: 4 : '((gr x) % %:o.x) % (1+x%~*:y)^-:>:x'

NB. (gr y) = (gamma 0.5*1+y) % gamma 0.5*y
gr=: 3 : 0
 select. (1=y)+2|y 
  case. 0 do. 0.5p0.5 * */%/ 0.5  0 +/ }.i.-:y    NB. even 
  case. 1 do. 2p_0.5  * */%/ 1  0.5 +/ }.i.-:y-1  NB. odd
  case. 2 do. 1p_0.5                              NB. 1=y  (% gamma 0.5)
 end.
)



----- Original Message -----
From: Roger Hui <rhui000 at shaw.ca>
Date: Thursday, August 14, 2008 11:21
Subject: Re: [Jprogramming] general Gamma distribution
To: Programming forum <programming at jsoftware.com>

> tpdf is the original formulation.  tpdf1 implements the 
> idea 
> in the "Alternatively" paragraph.
> 
> tpdf =: 4 : '((1+ x %~ *:y)^- -: >: x) % (0.5 beta -:x)*%: x'
> tpdf1=: 4 : '(gr x) % (%:o.x)*(1+x%~*:y)^-:>:x'
> 
> g15=: gamma 1.5
> 
> NB. (gr y) = (gamma 0.5*1+y) % gamma 0.5*y
> gr=: 3 : 0
>  select. (1=y)+2|y 
>   case. 0 do. g15 *  */%/ 0.5  0 +/ }.i.-
> :y    NB. even 
>   case. 1 do. g15 %~ */%/ 1  0.5 +/ }.i.-:y-1  
> NB. odd
>   case. 2 do. 
> 1p_0.5                           NB. 1=y  (% gamma 0.5)
>  end.
> )
> 
>    10 (tpdf -: tpdf1) 5 * 100 ?@$ 0
> 1
>    11 (tpdf -: tpdf1) 5 * 100 ?@$ 0
> 1
>    400 tpdf 1.6
> |NaN error: beta
> |   ((1+x%~*:y)^--:>:x)%(0.5     
> beta-:x)*%:x
>    400 tpdf1 1.6
> 0.11095
> 
> 
> 
> ----- Original Message -----
> From: Roger Hui <rhui000 at shaw.ca>
> Date: Thursday, August 14, 2008 7:55
> Subject: Re: [Jprogramming] general Gamma distribution
> To: Programming forum <programming at jsoftware.com>
> 
> > > John's version relies on ! which is limited to !170.  
> > When the mean exceeds
> > > 170, as Mr. Hui suggested,  perhaps the function should 
> > simply return a call
> > > to the normal function.
> > 
> > Alternatively, you could recode tpdf to use equation 5 in
> > http://mathworld.wolfram.com/Studentst-Distribution.html
> > and compute the ratio  (gamma 0.5*1+x)%gamma 0.5*x
> > not directly but by  */vec0%vec1 .
> > 
> > Separately, another possibility is to use the exact method 
> > up to its limit and apply integration on the remaining tail.
> > 
> > 
> > 
> > ----- Original Message -----
> > From: Robert Cyr <robert.cyr at gmail.com>
> > Date: Thursday, August 14, 2008 6:03
> > Subject: Re: [Jprogramming] general Gamma distribution
> > To: Programming forum <programming at jsoftware.com>
> > 
> > > Comparing the results of John's technique to Mr. Hui's exact 
> > > formula, it's
> > > interesting to note that
> > > 1. Simpson's rule produces results accurate to 2e_12 when 
> > > compared to the
> > > exact formulas.
> > > 2. the exact method is at least 3 times faster
> > > 3. they both have a problem with large mean
> > > 4. John's technique easily handles large x values
> > > 
> > > John's version relies on ! which is limited to !170.  
> > When 
> > > the mean exceeds
> > > 170, as Mr. Hui suggested,  perhaps the function should 
> > > simply return a call
> > > to the normal function.
> > > 
> > > Of course the last line of the int function refers to dp and 
> > not mp.
> > > 
> > > 
> > > On Wed, Aug 13, 2008 at 2:57 PM, John Randall 
> > > <randall at andromeda.rutgers.edu> wrote:
> > > 
> > > > Here is an alternative approach: numerically integrating 
> the 
> > > pdf.  It is
> > > > not subject to the bound in the beta function, and gives 
> > > sensible values
> > > > for large x and y.  The downside is figuring how to 
> get 
> > > the required
> > > > accuracy.  For the normal distribution, you can get 4-
> > > digit accuracy by
> > > > using h=0.05: I have used h=0.01.
> > > >
> > > > gamma=:!@:<:
> > > > beta =:gamma@[ * gamma@] % gamma at +
> > > >
> > > > NB. Probability density function
> > > > tpdf=:4 : '((1+ x %~ *:y)^- -: >: x) % (0.5 beta -:x)*%: x'
> > > >
> > > > ppr=: +//.@(*/)  NB. polynomial multiplication
> > > > dp=:+/ .*        NB. 
> dot 
> > product> >
> > > > NB. Simson's rule integrator
> > > > int=:1 : 0
> > > > n=.2>.+:>.y%0.01
> > > > h=.y%n
> > > > c=.1 4 1 ppr 0=2 | i. <: n
> > > > (h%3)*c mp u h*i.>:n
> > > > )
> > > >
> > > > NB. Find cdf by numerically integrating pdf
> > > > tcdf2=:4 : 0
> > > > if. y<0 do. 1-x tcdf -y return. end.
> > > > if. y=0 do. 0.5 return. end.
> > > > 0.5+x&tpdf int y
> > > > )
> > > >
> > > >
> > > >   5 tcdf2 1.47588
> > > > 0.899999
> > > >   5 tcdf2 2.01505
> > > > 0.95
> > > >   100 tcdf2 9
> > > > 1