[Rd] unstable corner of parameter space for qbeta?
peter dalgaard
pd@|gd @end|ng |rom gm@||@com
Thu Mar 26 09:09:57 CET 2020
It's a pretty extreme case, try e.g. curve(pbeta(x, shape1, shape2), n=10001), and (probably -- I can't be bothered to work out the relation between beta shapes and F df parameters this morning...) outside what is normally encountered in statistical analyses. Notice though, that you have lower=FALSE in your Qbeta0, so you should have it in qbeta as well:
> qbeta(t,shape1,shape2, lower=FALSE)
[1] 0.9999949
Warning message:
In qbeta(t, shape1, shape2, lower = FALSE) :
full precision may not have been achieved in 'qbeta'
which of course is still wrong (I dont't think there is a problem in the other tail, Qbeta0(t,...., lower=TRUE) returns zero.
You can see the effect also with
curve(qbeta(x, shape1, shape2), n=10001, from=.99, to=1)
which kind of suggests one of those regime-switching bugs, where different methods are used for various parts of the domain, and the cross-over is not done quite right.
At any rate, qbeta is one of R's very basic workhorses, so we do want it to work right, so by all means file a report.
-pd
> On 26 Mar 2020, at 02:09 , Ben Bolker <bbolker using gmail.com> wrote:
>
>
> I've discovered an infelicity (I guess) in qbeta(): it's not a bug,
> since there's a clear warning about lack of convergence of the numerical
> algorithm ("full precision may not have been achieved"). I can work
> around this, but I'm curious why it happens and whether there's a better
> workaround -- it doesn't seem to be in a particularly extreme corner of
> parameter space. It happens, e.g., for these parameters:
>
> phi <- 1.1
> i <- 0.01
> t <- 0.001
> shape1 = i/phi ## 0.009090909
> shape2 = (1-i)/phi ## 0.9
> qbeta(t,shape1,shape2) ## 5.562685e-309
> ## brute-force uniroot() version, see below
> Qbeta0(t,shape1,shape2) ## 0.9262824
>
> The qbeta code is pretty scary to read: the warning "full precision
> may not have been achieved" is triggered here:
>
> https://github.com/wch/r-source/blob/f8d4d7d48051860cc695b99db9be9cf439aee743/src/nmath/qbeta.c#L530
>
> Any thoughts? Should I report this on the bug list?
>
>
> A more general illustration:
> http://www.math.mcmaster.ca/bolker/misc/qbeta.png
>
> ===
> fun <- function(phi,i=0.01,t=0.001, f=qbeta) {
> f(t,shape1=i/phi,shape2=(1-i)/phi, lower.tail=FALSE)
> }
> ## brute-force beta quantile function
> Qbeta0 <- function(t,shape1,shape2,lower.tail=FALSE) {
> fn <- function(x) {pbeta(x,shape1,shape2,lower.tail=lower.tail)-t}
> uniroot(fn,interval=c(0,1))$root
> }
> Qbeta <- Vectorize(Qbeta0,c("t","shape1","shape2"))
> curve(fun,from=1,to=4)
> curve(fun(x,f=Qbeta),add=TRUE,col=2)
>
> ______________________________________________
> R-devel using r-project.org mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel
--
Peter Dalgaard, Professor,
Center for Statistics, Copenhagen Business School
Solbjerg Plads 3, 2000 Frederiksberg, Denmark
Phone: (+45)38153501
Office: A 4.23
Email: pd.mes using cbs.dk Priv: PDalgd using gmail.com
More information about the R-devel
mailing list