There appears to be a sign problem in the formula. The formula either form Wolfram or wikipedia are:
tcdf=: 4 : 0 assert. (%:x)>|y 0.5 + y * (!-:x+1) * ((0.5,-:1+x) H. 1.5 x%~-*:y) % (%:o.x) * !<:-:x ) As noted by others, the limit also needs to be adjusted. The limits probably are too high, as tcdf cannot exceed 1. It seems to be in that region that problem results occur. On Tue, Aug 12, 2008 at 5:38 PM, Roger Hui <[EMAIL PROTECTED]> wrote: > For large x one may want to use the N(0,1) CDF > in a pinch. From > http://www.jsoftware.com/jwiki/Essays/Normal_CDF > > erf =: (1 H. 1.5)@*: * 2p_0.5&* % ^@:*: > n01cdf=: -: @ >: @ erf @ %&(%:2) > > 10 tcdf 3 > 0.993328 > n01cdf 3 > 0.99865 > > The bad 100 tcdf 9 result is likely due to extreme > accumulation of numerical errors. On my machine I get: > > 100 tcdf 9 > 1.03437e17 > > n01cdf 9 > 1 > 1 - n01cdf 9 > _2.22045e_16 > > > > ----- Original Message ----- > From: John Randall <[EMAIL PROTECTED]> > Date: Tuesday, August 12, 2008 14:07 > Subject: Re: [Jprogramming] general Gamma distribution > To: Programming forum <[email protected]> > > > Roger Hui wrote: > > > It's the absolute value of y that should be bounded, > > > not just y itself. Thus: > > > > > > tcdf=: 4 : 0 > > > assert. (%:x)>|y > > > 0.5 + y * (!-:x-1) * ((0.5,-:1+x) H. 1.5 x%~-*:y) % > > (%:o.x) * !<:-:x > > > ) > > > > > > > > Mathworld <http://mathworld.wolfram.com/Studentst- > > Distribution.html> has > > different formulations (equations 7 and 8), which rely on > > the incomplete > > beta function and are not subject to the hypergeometric bounds. > > > > I don't think the restriction is bad when x (degrees of freedom) > > is above > > about 10. For example > > > > 10 tcdf 3 > > 0.993328 > > > > but it is not good if x is smaller. High values of y > > give strange results: > > > > 100 tcdf 9 > > _7.15204e15 > ---------------------------------------------------------------------- > For information about J forums see http://www.jsoftware.com/forums.htm > ---------------------------------------------------------------------- For information about J forums see http://www.jsoftware.com/forums.htm
