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.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.
