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