>Date: Fri, 2 May 2003 10:03:23 -0400 (EDT)
From: Morten Welinder <welinder@rentec.com>>To: p.dalgaard@biostat.ku.dk
>CC: r-devel@stat.math.ethz.ch, R-bugs@biostat.ku.dk
>Subject: Re: [Rd] qbeta hang (PR#2894)
>
>Ok, I can confirm that it does not, in fact, loop forever. Just a close
>approximation.
...
>There are lots of other places that worry me with respect to cancellation
>errors, for example
>
> r = 1 - pp;
> t = 1 - qq;
These can also be removed by changing qbeta.c:156-157 (R-1.7.1)
to
y = (y-a) * xinbta * (1.0-xinbta) *
exp (logbeta - pp * log(xinbta) - qq * log1p(-xinbta));
I don't understand why you have qbeta.c:167
if (fabs(y)<=acu) goto L_converged
which will fail for certain values of p & q, when x is close to zero as
the beta density tends to infinity. Why not use an exit condition based on
how far away from 'a' you are?
qbeta.c:174 (prevent excess precision)
Might this not just get optimized out? I doubt modern compilers respect
the volatile keyword. You should probably replace this with a safe
comparison. if(fabs(tx-xinbta)<=DBL_EPSILON*fabs(xinbta))...
** This should be checked with someone familiar with floating point
arithmetic; I'm not and may have got the correct expression wrong.
The 'infinite loop' is due to the speed of the pbeta routine, rather
than
the initial approximation from the qbeta. There is another algorithm
available for pbeta (Algorithm 708, CACM 18. 360--373, 1992) which may be
of use here. This continued fraction approximation is recommended by
Numerical Recipes over the power series approximation (make of that what
you will!)
TimM.