[Rd] unstable corner of parameter space for qbeta?

[peter dalgaard]

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