From jerry.lewis at biogenidec.com  Sun Jan  3 22:15:12 2010
From: jerry.lewis at biogenidec.com (jerry.lewis at biogenidec.com)
Date: Sun,  3 Jan 2010 22:15:12 +0100 (CET)
Subject: [Rd] dchisq tail accuracy (PR#14105)
Message-ID: <20100103211512.A79BA282EFCC@mail.pubhealth.ku.dk>

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]]



[Rd] dchisq tail accuracy (PR#14105)
jerry.lewis at biogenidec.com