unstable corner of parameter space for qbeta?

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

unstable corner of parameter space for qbeta?

Thread (9 messages)