On page 132 of the classic "Sharp APL Reference Manual"
by Paul Berry:
The APL function x!y is related to the complete Beta
function by the following identity:
Beta(x,y) <-> %y*(x-1)!x+y-1
Therefore, John Randall's original tpdf function can be
extended regarding its usefulness by:
beta=: ] [EMAIL PROTECTED] [ !&<: +
tpdf =: 4 : '((1+ x %~ *:y)^- -: >: x) % (0.5 beta -:x)*%: x'
11 (tpdf -: tpdf2) 10 * _0.5+100 [EMAIL PROTECTED] 0
1
10 (tpdf -: tpdf2) 10 * _0.5+100 [EMAIL PROTECTED] 0
1
400 (tpdf -: tpdf2) 1.6
1
----- Original Message -----
From: Roger Hui <[EMAIL PROTECTED]>
Date: Thursday, August 14, 2008 11:42
Subject: Re: [Jprogramming] general Gamma distribution
To: Programming forum <[email protected]>
> 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 <[EMAIL PROTECTED]>
> Date: Thursday, August 14, 2008 11:21
> Subject: Re: [Jprogramming] general Gamma distribution
> To: Programming forum <[email protected]>
>
> > 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 [EMAIL PROTECTED] 0
> > 1
> > 11 (tpdf -: tpdf1) 5 * 100 [EMAIL PROTECTED] 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 <[EMAIL PROTECTED]>
> > Date: Thursday, August 14, 2008 7:55
> > Subject: Re: [Jprogramming] general Gamma distribution
> > To: Programming forum <[email protected]>
> >
> > > > 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 <[EMAIL PROTECTED]>
> > > Date: Thursday, August 14, 2008 6:03
> > > Subject: Re: [Jprogramming] general Gamma distribution
> > > To: Programming forum <[email protected]>
> > >
> > > > 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
> > > > <[EMAIL PROTECTED]> 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 =:[EMAIL PROTECTED] * [EMAIL PROTECTED] % [EMAIL PROTECTED]
> > > > >
> > > > > 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
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
