You are dealing with functions that are non-zero over a very small interval, so you have to be very careful. There is no method that is going to be totally foolproof. Having said that, I have always felt that the default tolerance in integrate is too liberal (i.e. too large). I always use rel.tol of 1.e-08 (roughly, sqrt(machine epsilon)) in my computations, and I also increase subdivisions to 500.
Ravi. From: baptiste Auguié [mailto:baptiste.aug...@googlemail.com] Sent: Tuesday, September 21, 2010 9:58 AM To: Ravi Varadhan Cc: 'baptiste auguie'; 'r-help' Subject: Re: [R] puzzle with integrate over infinite range Hi, Thanks for the tip, but it's still mysterious to me. Reading ?integrate did not give me much hint as to what "relative accuracy" means in this context. I looked at the source of integrate.c but it's still not clear to me how the default value of rel.tol (10^-4 for me) is not enough to prevent a completely wrong answer (the error is much larger than this). Obviously, I'm worried now that I may not always choose a good value of ref.tol, if picked arbitrarily without my understanding of what it means. Thanks, baptiste On Sep 21, 2010, at 3:38 PM, Ravi Varadhan wrote: 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-boun...@r-project.org [mailto:r-help-boun...@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@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. [[alternative HTML version deleted]]
______________________________________________ R-help@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.
