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@r-project.org mailing list https://stat.ethz.ch/mailman/listinfo/r-devel