From: <jerry.lewis_at_biogenidec.com>

Date: Sun, 03 Jan 2010 22:15:12 +0100 (CET)

} while (q >= 1 || term * q > (1-q)*eps || term > termlimit);

/* lower tail */

term = mid; df = dfmid; i = imax;

while ( i ){

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

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.
*