[Rd] error in qbeta (PR#1201)

[Peter Dalgaard BSA]

z.yang@ucl.ac.uk writes:

> Full_Name: Ziheng Yang
> Version: 1.3.1
> OS: Windows 98
> Submission from: (NULL) (172.136.54.89)
> 
> 
> I noticed that qbeta is sometimes wrong and the error is not even due to the
> beta parameters being too extreme.  I am calculating the quantiles corresponding
> to cdf = 0.05, 0.15, ..., 0.95.  The value corresponding to cdf=0.25 is wrong
> while all other values are correct.  
> 
> qbeta(0.05, 0.143891, 0.05) = 1.040850e-05   correct
> qbeta(0.15, 0.143891, 0.05) = 0.021161       correct
> qbeta(0.25, 0.143891, 0.05) = 3e-308         wrong (correct value is 0.457227)
> qbeta(0.25, 0.143891, 0.05) = 0.945217       correct
> 
> I also used your source file qbeta.c to do the same calculation, and it seems
> that something might be wrong in the "initial approximation", before the
> newton-raphson iteration.  Also I notice that the 5/6 in the following line is
> probably supposed to mean 5./6., although this does not seem to be the source of
> the problem.
> 
> 	w = y * sqrt(h + r) / h - (t - s) * (r + 5 / 6 - 2 / (3 * h));
> 
> Could you possibly let me know how to fix the problem.  I will appreciate any
> help you might provide.  I am using the Windows version 1.3.1 (2001).

The 5/6 thing appears to be fixed in the current sources. However, the 
qbeta(0.25, 0.143891, 0.05) effect is still there and seems to be
algorithmic. It does not disappear when removing optimizations.

What happens is that in the calculation of the initial approximation,
we have 

        t = r * pow(1 - t + y * sqrt(t), 3);
        if (t <= 0.)
            xinbta = 1 - exp((log((1 - a) * qq) + logbeta) / qq);
... 

with t negative at the if() *and* xinbta coming out negative
so xinbta gets set to "lower" aka 3e-308

Further down, we get

   y = pbeta_raw(xinbta, pp, qq, /*lower_tail = */ TRUE);

which is way off from the target a, but then we calculate

   y = (y - a) *
         exp(logbeta + r * log(xinbta) + t * log(1 - xinbta));

which is a small number since log(xinbta) is large and negative. Then,
we have
 
   if (fabs(y) <= acu) goto L_converged;

and the whole thing terminates without even trying to get an improved
estimate.... 

Clues on where to fix this would be most welcome!

-- 
   O__  ---- Peter Dalgaard             Blegdamsvej 3  
  c/ /'_ --- Dept. of Biostatistics     2200 Cph. N   
 (*) \(*) -- University of Copenhagen   Denmark      Ph: (+45) 35327918
~~~~~~~~~~ - (p.dalgaard@biostat.ku.dk)             FAX: (+45) 35327907
-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-devel mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-devel-request@stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._