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