> 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 > > > > Best wishes, > > > > John ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
