[Rd] dchisq tail accuracy (PR#14105)
jerry.lewis at biogenidec.com
jerry.lewis at biogenidec.com
Sun Jan 3 22:15:12 CET 2010
The issue seems to be that the infinite sum is truncated too early when x
is in the extreme upper tail. An easily validated improvement to to
dnchisq.c would be to add an additional requirement in the upper tail
while condition, that the summation should continue while the additive
term remains too big a fraction of the saddlepoint density. Here is
replacement code that partially offsets the time required for addtional
terms in the sum by calculating x*ncp2 only once instead of repeating it
in every iteration of both loops.
x2 = x * ncp2
/* upper tail */
term = mid; df = dfmid; i = imax;
s = 0.5 - (dx+sqrt(dx*dx+4*ncp/x))*0.25 /* saddlepoint */
s2 = 1-2*s
termlimit = exp(ncp*s/s2 - log(s2)*df*0.5 - s*x -
log((df+2*ncp+4*ncp*s/s2)/(s2*s2)*pi) - log(2)*2**54) /* saddlepoint
density * 2^-53 */
do {
i++;
q = x2 / i / df;
df += 2;
term *= q;
sum += term;
} while (q >= 1 || term * q > (1-q)*eps || term > termlimit);
/* lower tail */
term = mid; df = dfmid; i = imax;
while ( i ){
df -= 2;
q = i * df / x2;
i--;
term *= q;
sum += term;
if (q < 1 && term * q <= (1-q)*eps) break;
}
Jerry
[[alternative HTML version deleted]]
More information about the R-devel
mailing list