

I have been told by the CRAN administrators that the following code generated
an error on 64bit 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.5e8
> 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 1dimensional
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 pre2.12.0 behaviour. (PR#15219)",
then I do not understand what this pre2.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/rdevel


On 20130716 07:55, Hans W Borchers wrote:
> I have been told by the CRAN administrators that the following code generated
> an error on 64bit 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.5e8
> > 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 1dimensional
> 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 pre2.12.0 behaviour. (PR#15219)",
>
> then I do not understand what this pre2.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.5e8
> 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.233257e07
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/rdevel


On Tue, 20130716 at 13:55 +0200, Hans W Borchers wrote:
> I have been told by the CRAN administrators that the following code generated
> an error on 64bit 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.5e8
> > 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 1dimensional
> 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 pre2.12.0 behaviour. (PR#15219)",
>
> then I do not understand what this pre2.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.rproject.org/bugzilla3/show_bug.cgi?id=15219It 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 misreported.
The tolerance in your example (1.5e8) is considerably smaller than the
default (1.2e4). 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.5e7 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/rdevel
This message and its attachments are strictly confidenti...{{dropped:8}}
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rdevel


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, 20130716 at 13:55 +0200, Hans W Borchers wrote:
>> I have been told by the CRAN administrators that the following code generated
>> an error on 64bit 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.5e8
>> > 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 1dimensional
>> 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 pre2.12.0 behaviour. (PR#15219)",
>>
>> then I do not understand what this pre2.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.rproject.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 misreported.
>
> The tolerance in your example (1.5e8) is considerably smaller than the
> default (1.2e4). 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.5e7 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/rdevel>
> 
> This message and its attachments are strictly confiden...{{dropped:8}}
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rdevel


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, 20130716 at 13:55 +0200, Hans W Borchers wrote:
>> I have been told by the CRAN administrators that the following code
>> generated an error on 64bit 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.5e8
>> > 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
>> 1dimensional 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 pre2.12.0 behaviour. (PR#15219)",
>>
>> then I do not understand what this pre2.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.rproject.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 misreported.
>
> The tolerance in your example (1.5e8) is considerably smaller than
> the default (1.2e4). 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.5e7 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/rdevel>
> 
>  This message and its attachments are strictly
> confiden...{{dropped:8}}
______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rdevel______________________________________________
[hidden email] mailing list
https://stat.ethz.ch/mailman/listinfo/rdevel

