Problem following an R bug fix to integrate()

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 small
error 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}}

Problem following an R bug fix to integrate()

Thread (5 messages)