>Date: Fri, 2 May 2003 10:03:23 -0400 (EDT) From: Morten Welinder <[EMAIL PROTECTED]> >To: [EMAIL PROTECTED] >CC: [EMAIL PROTECTED], [EMAIL PROTECTED] >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. ______________________________________________ [EMAIL PROTECTED] mailing list https://www.stat.math.ethz.ch/mailman/listinfo/r-devel