[Rd] unstable corner of parameter space for qbeta?

[Ben Bolker]

On 2020-03-26 4:02 a.m., Martin Maechler wrote:
>>>>>> Ben Bolker 
>>>>>>     on Wed, 25 Mar 2020 21:09:16 -0400 writes:
> 
>     > 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?
> 
> Well,  qbeta() is mostly based on inverting pbeta()  and pbeta()
> has *several* "dangerous" corners in its parameter spaces
> {in some cases, it makes sense to look at the 4 different cases
>  log.p = TRUE/FALSE  //  lower.tail = TRUE/FALSE  separately ..}
> 
> pbeta() itself is based on the most complex numerical code in
> all of base R, i.e., src/nmath/toms708.c  and that algorithm
> (TOMS 708) had been sophisticated already when it was published,
> and it has been improved and tweaked several times since being
> part of R, notably for the log.p=TRUE case which had not been in
> the focus of the publication and its algorithm.
> [[ NB: part of this you can read when reading  help(pbeta)  to the end ! ]]
> 
> I've spent many "man weeks", or even "man months" on pbeta() and
> qbeta(), already and have dreamed to get a good student do a
> master's thesis about the problem and potential solutions I've
> looked into in the mean time.
> 
> My current gut feeling is that in some cases, new approximations
> are necessary (i.e. tweaking of current approximations is not
> going to help sufficiently).
> 
> Also not (in the R sources)  tests/p-qbeta-strict-tst.R
> a whole file of "regression tests" about  pbeta() and qbeta()
> {where part of the true values have been computed with my CRAN
> package Rmpfr (for high precision computation) with the
> Rmpfr::pbetaI() function which gives arbitrarily precise pbeta()
> values but only when  (a,b) are integers -- that's the "I" in pbetaI().
> 
> Yes, it's intriguing ... and I'll look into your special
> findings a bit later today.
> 
> 
>   > Should I report this on the bug list?
> 
> Yes, please.  Not all problem of pbeta() / qbeta() are part yet,
> of R's bugzilla data base,  and maybe this will help to draw
> more good applied mathematicians look into it.

  Will report.

  I'm not at all surprised that this is a super-tough problem.  The only
part that was surprising to me was that my naive uniroot-based solution
worked (for this particular corner of parameter space where qbeta() has
trouble: it was terrible elsewhere, so now I'm using a hybrid solution
where I use my brute-force uniroot thing if I get a warning from qbeta().

I hesitated to even bring it up because I know you're really busy, but I
figured it was better to tag it now and let you deal with it some time
later.

Bugzilla report at https://bugs.r-project.org/bugzilla/show_bug.cgi?id=17746

  cheers
   Ben Bolker


> 
> 
> 
> Martin Maechler
> ETH Zurich and R Core team
> (I'd call myself the "dpq-hacker" within R core -- related to
>  my CRAN package 'DPQ')
> 
> 
>     > 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
>