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.
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
On Tue, 2013-07-16 at 13:55 +0200, 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.You can see the bug report here: https://bugs.r-project.org/bugzilla3/show_bug.cgi?id=15219 It concerns the behaviour of integrate with a small error tolerance.>From 2.12.0 to 3.0.1 integrate was not working correctly with smallerror tolerance values, in the sense that small values did not improve accuracy and the accuracy was mis-reported. The tolerance in your example (1.5e-8) is considerably smaller than the default (1.2e-4). My guess is that the rounding error always existed but was not detected due to the bug. You might try a larger tolerance. I have tested your example and increasing the tolerance to 1.5e-7 removes the error. Martyn> 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. > > ______________________________________________ > R-devel at r-project.org mailing list > https://stat.ethz.ch/mailman/listinfo/r-devel----------------------------------------------------------------------- This message and its attachments are strictly confidenti...{{dropped:8}}