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
