Problem following an R bug fix to integrate()

classic Classic list List threaded Threaded
5 messages Options
Reply | Threaded
Open this post in threaded view
|

Problem following an R bug fix to integrate()

Hans W Borchers-2
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.

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Problem following an R bug fix to integrate()

J. R. M. Hosking-3
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

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Problem following an R bug fix to integrate()

Martyn Plummer-3
In reply to this post by Hans W Borchers-2
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 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.
>
> ______________________________________________
> [hidden email] mailing list
> https://stat.ethz.ch/mailman/listinfo/r-devel

-----------------------------------------------------------------------
This message and its attachments are strictly confidenti...{{dropped:8}}

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Problem following an R bug fix to integrate()

Hans W Borchers-2
Thanks for the help.

What bothers me is that it works on most systems and does not work on
some more 'exotic' systems -- though it should work everywhere however
small the user chooses the tolerance (with some warnings, maybe).

I decided I will apply my own integration routines in this example as
they appear to work more reliably.

Hans Werner


On Wed, Jul 17, 2013 at 7:37 PM, Martyn Plummer <[hidden email]> wrote:

> 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 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.
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
> -----------------------------------------------------------------------
> This message and its attachments are strictly confiden...{{dropped:8}}

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel
Reply | Threaded
Open this post in threaded view
|

Re: Problem following an R bug fix to integrate()

Ravi Varadhan-2
This, i.e. quadrature, is another area where the "default" or "base" R functionality needs enhancement, just like the functionality for optimization.  While `integrate' is good, it can be improved.  

Hans Werner, what routines do you use for quadrature?

Best,
Ravi

-----Original Message-----
From: [hidden email] [mailto:[hidden email]] On Behalf Of Hans W Borchers
Sent: Wednesday, July 17, 2013 3:36 PM
To: Martyn Plummer
Cc: [hidden email]
Subject: Re: [Rd] Problem following an R bug fix to integrate()

Thanks for the help.

What bothers me is that it works on most systems and does not work on some more 'exotic' systems -- though it should work everywhere however small the user chooses the tolerance (with some warnings, maybe).

I decided I will apply my own integration routines in this example as they appear to work more reliably.

Hans Werner


On Wed, Jul 17, 2013 at 7:37 PM, Martyn Plummer <[hidden email]> wrote:

> 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 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.
>>
>> ______________________________________________
>> [hidden email] mailing list
>> https://stat.ethz.ch/mailman/listinfo/r-devel
>
> ----------------------------------------------------------------------
> - This message and its attachments are strictly
> confiden...{{dropped:8}}

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel

______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/r-devel