[Rd] Problem following an R bug fix to integrate()
J. R. M. Hosking
JRMHosking at gmail.com
Wed Jul 17 16:38:55 CEST 2013
On 2013-07-16 07:55, Hans W Borchers wrote:
> I have been told by the CRAN administrators that the following code generated
> an error on 64-bit Fedora Linux (gcc, clang) and on Solaris machines (sparc,
> x86), but runs well on all other systems):
>
> > fn<- function(x, y) ifelse(x^2 + y^2<= 1, 1 - x^2 - y^2, 0)
>
> > tol<- 1.5e-8
> > fy<- function(x) integrate(function(y) fn(x, y), 0, 1,
> subdivisions = 300, rel.tol = tol)$value
> > Fy<- Vectorize(fy)
>
> > xa<- -1; xb<- 1
> > Q<- integrate(Fy, xa, xb,
> subdivisions = 300, rel.tol = tol)$value
>
> Error in integrate(Fy, xa, xb, subdivisions = 300, rel.tol = tol) :
> roundoff error was detected
>
> Obviously, this realizes a double integration, split up into two 1-dimensional
> integrations, and the result shall be pi/4. I wonder what a 'roundoff error'
> means in this situation.
>
> In my package, this test worked well, w/o error or warnings, since July 2011,
> on Windows, Max OS X, and Ubuntu Linux. I have no chance to test it on one of
> the above mentioned systems. Of course, I can simply disable these tests, but
> I would not like to do so w/o good reason.
>
> If there is a connection to a bug fix to integrate(), with NEWS item
>
> "integrate() reverts to the pre-2.12.0 behaviour. (PR#15219)",
>
> then I do not understand what this pre-2.12.0 behavior really means.
>
> Thanks for any help or a hint to what shall be changed.
> Hans W Borchers
>
> PS:
> This kind of tricky definition in function 'fn' has caused some discussion on
> this list in July 2009. I still think it should be allowed to proceed in this
> way.
Short answer: use a larger value of 'rel.tol' for the outer integral
than for the inner integral.
Using R 2.11.1 on Windows:
> fn <- function(x, y) ifelse(x^2 + y^2 <= 1, 1 - x^2 - y^2, 0)
> tol <- 1.5e-8
> fy <- function(x) integrate(function(y) fn(x, y), 0, 1,
subdivisions = 300, rel.tol = tol)$value
> Fy <- Vectorize(fy)
> xa <- -1; xb <- 1
> Q <- integrate(Fy, xa, xb,
subdivisions = 300, rel.tol = tol)$value
Error in integrate(Fy, xa, xb, subdivisions = 300, rel.tol = tol) :
roundoff error was detected
Now increase 'rel.tol' in the outer integral:
> Q <- integrate(Fy, xa, xb,
subdivisions = 300, rel.tol = tol*10)$value
> Q - pi/4
[1] -1.233257e-07
Longer answer: Fy, the integrand of the outer integral, is in effect
computed with noise added to it that is of the order of magnitude of
the 'rel.tol' of the inner integral; this noise prevents the outer
integral from attaining relative accuracy of this magnitude or smaller.
The version of integrate() in use since R 2.12.0 did not accurately
reproduce the computations in the Fortran routines (in the QUADPACK
package) on which it was based, and in consequence failed to detect this
situation. Reversion to the R 2.11.1 version of integrate() restores
concordance with the Fortran routines and correctly diagnoses the
inability of the outer integral to achieve the requested accuracy.
(And, btw, the Q computed above is actually closer to pi/4 than you
will have been getting with the code that "worked well".)
J. R. M. Hosking
More information about the R-devel
mailing list