Re: [Rd] dchisq tail accuracy (PR#14105)

From: <jerry.lewis_at_biogenidec.com>
Date: Sun, 03 Jan 2010 22:15:12 +0100 (CET)


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



R-devel_at_r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel Received on Sun 03 Jan 2010 - 22:45:08 GMT

Archive maintained by Robert King, hosted by the discipline of statistics at the University of Newcastle, Australia.
Archive generated by hypermail 2.2.0, at Mon 04 Jan 2010 - 10:20:10 GMT.

Mailing list information is available at https://stat.ethz.ch/mailman/listinfo/r-devel. Please read the posting guide before posting to the list.

list of date sections of archive