Bug 15755 - qbeta wrong behavior
Summary: qbeta wrong behavior
Status: CLOSED FIXED
Alias: None
Product: R
Classification: Unclassified
Component: Accuracy (show other bugs)
Version: R 3.1.0
Hardware: x86_64/x64/amd64 (64-bit) Windows 64-bit
: P5 major
Assignee: Martin Maechler
URL:
Depends on:
Blocks:
 
Reported: 2014-04-17 12:05 UTC by edern le roux
Modified: 2015-12-14 13:45 UTC (History)
5 users (show)

See Also:


Attachments

Note You need to log in before you can comment on or make changes to this bug.
Description edern le roux 2014-04-17 12:05:49 UTC
With some values of p, shape1 and shape2, qbeta returns 1 while it should return something around 0:

> qbeta(0.6948886,0.0672788,226390)
[1] 1
Warning message:
In qbeta(0.6948886, 0.0672788, 226390) :
  full precision may not have been achieved in 'qbeta'
Comment 1 Roby Joehanes 2014-04-29 20:57:35 UTC
Could this be related to bug #15641? If so, the fix is in the next R release.
Comment 2 Roby Joehanes 2014-05-02 15:35:04 UTC
This bug apparently still affects the latest R-devel or R-patched. I also noticed that qbeta(0.6948886,226390, 0.0672788) == 1, which is correct (true answer is 1-7.884458419482598E-9). So, I guess the solution is to swap p and q when needed (in qbeta.c).

Here is my proposed patch in qbeta.c:

    /* calculate the initial approximation */
if (q / p > 1e7) {
double temp = pp; pp = qq; qq = temp;
swap_tail = !swap_tail;
}
    /* y := {fast approximation of} qnorm(1 - a) :*/
Comment 3 Martin Maechler 2014-05-02 16:14:35 UTC
(In reply to Roby Joehanes from comment #2)
> This bug apparently still affects the latest R-devel or R-patched. 

yes, I know.
I don't want to adopt another fudge factor ("1e7") solution.
Rather, I'm aiming for switching to use the Newton approximation to
"log log scale" in such extreme cases.  This will be much faster converging
and I think not needing the inner loop for Newton step size adjustment
(which has been slightly dubious in my view anyway, also the test that Newton step sizes must montonely decrease etc, seems adhoc and only needed because in such situations, "original scale Newton" seems inappropriate.
Comment 4 Umanga 2014-10-08 09:55:57 UTC
(In reply to edern le roux from comment #0)
> With some values of p, shape1 and shape2, qbeta returns 1 while it should
> return something around 0:
> 
> > qbeta(0.6948886,0.0672788,226390)
> [1] 1
> Warning message:
> In qbeta(0.6948886, 0.0672788, 226390) :
>   full precision may not have been achieved in 'qbeta'

Testttttttttttttttttttttttttttttthjhdshfhdshfdsfhdsdsfdfhjshjghghdfhgdf
Comment 5 Martin Maechler 2015-02-26 15:44:19 UTC
(In reply to Martin Maechler from comment #3)
> (In reply to Roby Joehanes from comment #2)
> > This bug apparently still affects the latest R-devel or R-patched. 
> 
> yes, I know.
> I don't want to adopt another fudge factor ("1e7") solution.
> Rather, I'm aiming for switching to use the Newton approximation to
> "log log scale" in such extreme cases.  This will be much faster converging
> and I think not needing the inner loop for Newton step size adjustment
> (which has been slightly dubious in my view anyway, also the test that
> Newton step sizes must montonely decrease etc, seems adhoc and only needed
> because in such situations, "original scale Newton" seems inappropriate.

I have finally committed (svn rev 67902) a much extended 
  src/nmath/qbeta.c
which fixes this bug and quite a few others.
Indeed, now using  Newton on log scale  or log-log scale,
and doing quite a few other careful checks etc.

The new code is not at all beautiful and definitely more baroque than the previous code.  OTOH, these changes have been waiting now for many months, and have never failed in my usage and tests.

So they should become part of R 3.2.x (and later)... and may prompt more beautiful code submissions which works in as many cases as this one.
Comment 6 Ben Kearns 2015-03-02 13:32:25 UTC
(In reply to Martin Maechler from comment #5)
> (In reply to Martin Maechler from comment #3)
> > (In reply to Roby Joehanes from comment #2)
> > > This bug apparently still affects the latest R-devel or R-patched. 
> > 
> > yes, I know.
> > I don't want to adopt another fudge factor ("1e7") solution.
> > Rather, I'm aiming for switching to use the Newton approximation to
> > "log log scale" in such extreme cases.  This will be much faster converging
> > and I think not needing the inner loop for Newton step size adjustment
> > (which has been slightly dubious in my view anyway, also the test that
> > Newton step sizes must montonely decrease etc, seems adhoc and only needed
> > because in such situations, "original scale Newton" seems inappropriate.
> 
> I have finally committed (svn rev 67902) a much extended 
>   src/nmath/qbeta.c
> which fixes this bug and quite a few others.
> Indeed, now using  Newton on log scale  or log-log scale,
> and doing quite a few other careful checks etc.
> 
> The new code is not at all beautiful and definitely more baroque than the
> previous code.  OTOH, these changes have been waiting now for many months,
> and have never failed in my usage and tests.
> 
> So they should become part of R 3.2.x (and later)... and may prompt more
> beautiful code submissions which works in as many cases as this one.

Glad to hear that this has been fixed, and will be shortly made available.

As a relative newbie to R, I was wondering if there's a simple way of incorporating the updated qbeta function earlier? For example, if the updated code's available code this be incorporated as a user-defined function for now?

Thanks
Comment 7 Martin Maechler 2015-03-17 16:05:21 UTC
(In reply to Ben Kearns from comment #6)
> (In reply to Martin Maechler from comment #5)
   [.................]
 
> > So they should become part of R 3.2.x (and later)... and may prompt more
> > beautiful code submissions which works in as many cases as this one.
> 
> Glad to hear that this has been fixed, and will be shortly made available.
> 
> As a relative newbie to R, I was wondering if there's a simple way of
> incorporating the updated qbeta function earlier? For example, if the
> updated code's available code this be incorporated as a user-defined
> function for now?

As always: You download an "R-devel snapshot" - in your case, for Windows,
from your favorite CRAN mirror, e.g., from
http://cran.rstudio.org/bin/windows/base/rdevel.html