[Rd] Integrate erros on certain functions

[John P. Nolan]

Dear John,

Two issues.  First, the default action is stop.on.error=TRUE, so anytime the integrate function can determine an error, it will stop.  It doesn't detect an error, so no error is produced (whether you set stop.on.error=TRUE or FALSE).

The second issue is the real problem:  what automatic numerical integration can do or not do.  When the upper bound is Inf, integrate (which is based on the QUADPACK fortran code) does a change of variable to make the region of integration be a finite interval, then tries to evaluate that transformed integral.  When the upper bound is a large finite number, integrate tries to evaluate the integral directly.  In this case, the integrand is evaluated at multiple x values in the interval [0,13000].   Those x values are large, and the resulting function values are very near 0.  You can verify this by putting some trace statements into your integrand function, e.g. 

> f1 <- function( x ) { y <- exp(-x); print( rbind(x,y) ); return(y) }
> integrate( f1, lower = 0, upper =13000)

The quadrature rule "sees" an integrand near 0 and returns a value for the integral near 0.  It does not detect an error, so it does not report anything to you.  It does not know how the integrand function behaves on regions where it does not evaluate it.  This is a well-known problem in numerical integration: there is no way the integrate function can know what region to focus on in a general problem.  Using an upper bound=Inf does not guarantee that you will get the correct value, but sometimes it works.  

Hope this helps. 

John

………………………………………………………………………………..
John P. Nolan
Math/Stat Dept., American University
106J Myers Hall, 4400 Massachusetts Ave, NW, Washington, DC 20016-8050
Phone: 202-885-3140   E-mail:  jpnolan at american.edu
Web:   http://fs2.american.edu/jpnolan/www/





-----Original Message-----
From: R-devel <r-devel-bounces at r-project.org> On Behalf Of John Muschelli
Sent: Friday, March 23, 2018 6:52 PM
To: r-devel at r-project.org
Subject: [Rd] Integrate erros on certain functions

In the help for ?integrate:

>When integrating over infinite intervals do so explicitly, rather than
just using a large number as the endpoint. This increases the chance of a correct answer – any function whose integral over an infinite interval is finite must be near zero for most of that interval.

I understand that and there are examples such as:

## a slowly-convergent integral
integrand <- function(x) {1/((x+1)*sqrt(x))} integrate(integrand, lower = 0, upper = Inf)

## don't do this if you really want the integral from 0 to Inf integrate(integrand, lower = 0, upper = 1000000, stop.on.error = FALSE) #> failed with message ‘the integral is probably divergent’

which gives an error message if stop.on.error = FALSE. But what happens on something like the function below:
integrate(function(x) exp(-x), lower = 0, upper =Inf) #> 1 with absolute error < 5.7e-05
integrate(function(x) exp(-x), lower = 0, upper =13000) #> 2.819306e-05 with absolute error < 5.6e-05

*integrate(function(x) exp(-x), lower = 0, upper =13000, stop.on.error = FALSE)#> 2.819306e-05 with absolute error < 5.6e-05*

I'm not sure this is a bug or misuse of the function, but I would assume the last integrate to give an error if stop.on.error = FALSE.

	[[alternative HTML version deleted]]

______________________________________________
R-devel at r-project.org mailing list
https://urldefense.proofpoint.com/v2/url?u=https-3A__stat.ethz.ch_mailman_listinfo_r-2Ddevel&d=DwIFaQ&c=U0G0XJAMhEk_X0GAGzCL7Q&r=7rQvU8hscCTWlvO-F5wI2-2eTiW40XI5qUKda0AnbG0&m=iA2KskSHO_cMznVT31Amx5mIJ0-cQurEM9ItQz-WwvU&s=_A2zZDw5gLetKaZqbPZMpJFqO8B1-kPT2T__T73CM-I&e=