[R] puzzle with integrate over infinite range

[Ravi Varadhan]

There is nothing mysterious.  You need to increase the accuracy of
quadrature by decreasing the error tolerance:

# I scaled your function to a proper Gaussian density
shiftedGauss <- function(x0=500){
 integrate(function(x) 1/sqrt(2*pi * 100^2) * exp(-(x-x0)^2/(2*100^2)), 0,
Inf, rel.tol=1.e-07)$value }

shift <- seq(500, 800, by=10)
plot(shift, sapply(shift, shiftedGauss))


Hope this helps,
Ravi.

-----Original Message-----
From: r-help-bounces at r-project.org [mailto:r-help-bounces at r-project.org] On
Behalf Of baptiste auguie
Sent: Tuesday, September 21, 2010 8:38 AM
To: r-help
Subject: [R] puzzle with integrate over infinite range

Dear list,

I'm calculating the integral of a Gaussian function from 0 to
infinity. I understand from ?integrate that it's usually better to
specify Inf explicitly as a limit rather than an arbitrary large
number, as in this case integrate() performs a trick to do the
integration better.

However, I do not understand the following, if I shift the Gauss
function by some amount the integral should not be affected,

shiftedGauss <- function(x0=500){
 integrate(function(x) exp(-(x-x0)^2/100^2), 0, Inf)$value
}

shift <- seq(500, 800, by=10)
plot(shift, sapply(shift, shiftedGauss))

Suddenly, just after 700, the value of the integral drops to nearly 0
when it should be constant all the way. Any clue as to what's going on
here? I guess it's suddenly missing the important part of the range
where the integrand is non-zero, but how could this be overcome?

Regards,

baptiste


sessionInfo()
R version 2.11.1 (2010-05-31)
x86_64-apple-darwin9.8.0

locale:
[1] en_GB.UTF-8/en_GB.UTF-8/C/C/en_GB.UTF-8/en_GB.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base

other attached packages:
[1] inline_0.3.5        RcppArmadillo_0.2.6 Rcpp_0.8.6
statmod_1.4.6

loaded via a namespace (and not attached):
[1] tools_2.11.1

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.