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
