[Rd] unstable corner of parameter space for qbeta?

Ravi Varadhan r@v|@v@r@dh@n @end|ng |rom jhu@edu
Thu Mar 26 19:33:43 CET 2020

This is also strange:

qbeta <- function (p, shape1, shape2, ncp = 0, lower.tail = TRUE, log.p = FALSE)
    if (missing(ncp))
        .Call(C_qbeta, p, shape1, shape2, lower.tail, log.p)
    else .Call(C_qnbeta, p, shape1, shape2, ncp, lower.tail,

Since the default value is 0 for non-centrality, it seems like the logic above is wrong. When ncp=0, C_qnbeta would be called rather than C_qbeta.

Am I missing something?


From: R-devel <r-devel-bounces using r-project.org> on behalf of J C Nash <profjcnash using gmail.com>
Sent: Thursday, March 26, 2020 10:40:05 AM
To: Martin Maechler
Cc: r-devel using r-project.org
Subject: Re: [Rd] unstable corner of parameter space for qbeta?

Despite the need to focus on pbeta, I'm still willing to put in some effort.
But I find it really helps to have 2-3 others involved, since the questions back
and forth keep matters moving forward. Volunteers?

Thanks to Martin for detailed comments.


On 2020-03-26 10:34 a.m., Martin Maechler wrote:
>>>>>> J C Nash
>>>>>>     on Thu, 26 Mar 2020 09:29:53 -0400 writes:
>     > Given that a number of us are housebound, it might be a good time to try to
>     > improve the approximation. It's not an area where I have much expertise, but in
>     > looking at the qbeta.c code I see a lot of root-finding, where I do have some
>     > background. However, I'm very reluctant to work alone on this, and will ask
>     > interested others to email off-list. If there are others, I'll report back.
> Hi John.
> Yes, qbeta() {in its "main branches"}  does zero finding, but
> zero finding of   pbeta(...) - p*   and I tried to explain in my
> last e-mail that the real problem is that already pbeta() is not
> accurate enough in some unstable corners ...
> The order fixing should typically be
> 1) fix pbeta()
> 2) look at qbeta() which now may not even need a fix because its
>    problems may have been entirely a consequence of pbeta()'s inaccuracies.
>    And if there are cases where the qbeta() problems are not
>    only pbeta's "fault", it is still true that the fixes that
>    would still be needed crucially depend on the detailed
>    working of the function whose zero(s) are sought, i.e.,  pbeta()
>     > Ben: Do you have an idea of parameter region where approximation is poor?
>     > I think that it would be smart to focus on that to start with.
> ----------------------------
> Rmpfr  matrix-/vector - products:
>     > Martin: On a separate precision matter, did you get my query early in year about double
>     > length accumulation of inner products of vectors in Rmpfr? R-help more or
>     > less implied that Rmpfr does NOT use extra length. I've been using David
>     > Smith's FM Fortran where the DOT_PRODUCT does use double length, but it
>     > would be nice to have that in R. My attempts to find "easy" workarounds have
>     > not been successful, but I'll admit that other things took precedence.
> Well, the current development version of 'Rmpfr' on R-forge now
> contains facilities to enlarge the precision of the computations
> by a factor 'fPrec' with default 'fPrec = 1';
> notably, instead of  x %*% y   (where the `%*%` cannot have more
> than two arguments) does have a counterpart  matmult(x,y, ....)
> which allows more arguments, namely 'fPrec', or directly 'precBits';
> and of course there are  crossprod() and tcrossprod() one should
> use when applicable and they also got the  'fPrec' and
> 'precBits' arguments.
> {The %*% etc precision increase still does not work optimally
>  efficiency wise, as it simply increases the precision of all
>  computations by just increasing the precision of x and y (the inputs)}.
> The whole  Matrix and Matrix-vector arithmetic is still
> comparibly slow in Rmpfr .. mostly because I valued human time
> (mine!) much higher than computer time in its implementation.
> That's one reason I would never want to double the precision
> everywhere as it decreases speed even more, and often times
> unnecessarily: doubling the accuracy is basically "worst-case
> scenario" precaution
> Martin

R-devel using r-project.org mailing list

	[[alternative HTML version deleted]]

More information about the R-devel mailing list