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