[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