Thanks, I'll do that too from now on. 
It strikes me that in a case such as this one it may be safer to use a truncated, finite interval around the region where the integrand is non-zero, rather than following the advice of ?integrate to use Inf as integration limit. At least one wouldn't risk to get an entirely wrong result depending on a choice of rel.tol. Regarding this parameter, is there a simple interpretation of how it affected the result in the context of my example?

Thanks again,

baptiste

On Sep 21, 2010, at 4:04 PM, Ravi Varadhan wrote:

> 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.auguie@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-bounces@r-project.org [mailto:r-help-bounces@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]]

