# Monte Carlo integration help needed

I'm trying to simulate these two integrals using Monte Carlo simulation: $$\int_{-\infty}^\infty \exp(-x^2) dx, \quad \mbox{and } \int_{-\infty}^\infty \exp(-|x|) dx .$$

When I use runif(n,-Inf,Inf), I get NaN (number very close to zero) (using R). I also tried converting that to a polar coordinates integral ($drdθ$ integral) but then I don't really know how to do the double integration needed.

Any ideas?

• The integral for $exp(-x^2)$ using polar coordinates is relatively easy if you work with the square of the integral $I^2 = \int exp(-x^2)dx.\int exp(-y^2)dy$; - write as a double integral, convert to polar (use $x^2+y^2=r^2$ along the way), don't forget the Jacobian, and you should be set. Commented Apr 15, 2014 at 1:40
• I already tried that, but I'm stuck at how to do a double integration in R, not on paper Commented Apr 15, 2014 at 12:49
• Oh. In that case (that you want to do the double integral numerically) I think the conversion to polar coordinates may not be the best choice (though the polar method for generating normal random variables is one of two methods that does use polar coordinates, and so could form the basis for such a method of working out that integral) Commented Apr 15, 2014 at 19:21

For suitably large $a>0$, you can sample from a $\mathrm{U}[-a,a]$ distribution because, for these two integrals, since the "tails" of the functions decay fast enough, $\int_{-\infty}^{-a} f(x)\,dx + \int_{a}^\infty f(x)\,dx$ is small. Don't forget to normalize the uniform density properly. For comparison, here are the analytical results and the corresponding simulations.

$$\int_{-\infty}^\infty e^{-x^2}\,dx = \sqrt{\pi} \approx 1.772454$$

N <- 10^6
a <- 5
(2*a) * sum(exp(-(runif(N, -a, a)^2))) / N
[1] 1.77192


$$\int_{-\infty}^\infty e^{-|x|}\,dx = 2$$

(2*a) * sum(exp(-abs(runif(N, -a, a)))) / N
[1] 1.988736

• So you use the uniform distribution but in the range of (-5,5), but still get the correct result ? Could you explain that please ? Commented Apr 15, 2014 at 12:13
• I see. But how come I get worse results if I even make that range (-6,6) ? The result diverges greatly. What makes (-5,5) the right one ? Also, is there a way to make a graph of this simulation (of what it does in every step, I guess)? Commented Apr 15, 2014 at 14:29
• It doesn't diverge: N <- 10^6; 12 * sum(exp(-(runif(N, -6, 6)^2))) / N gives 1.775803 for example. You probably didn't adjust the normalization constant of the density when you went from $[-5,5]$ to $[-6,6]$. That $10$ must become $12$.
– Zen
Commented Apr 15, 2014 at 16:25
• Although you are correct, your reasoning is not: the smallness of the integrand outside any finite interval does not guarantee that it can be ignored outside that interval!
– whuber
Commented Apr 15, 2014 at 17:09
• Awesome. It's just that as a mathematician, this (-5,5) seemed weird. I get it now with the alphas. Thanks a lot. Any idea about the graph ? I tried running the function for some N's with a for loop and a dotplot from lattice, but it didn't work. Commented Apr 15, 2014 at 18:45

According with runif manual page which you can find here, the min and max parameter values must be finite. If not it will produce NAs.

• This is an example of a poorly written code. It should have raised an exception or better a fatal error when called with Inf argument instead of silently returning NaN. Commented Apr 14, 2014 at 15:28
• @Aksakal Have you made this suggestion to the R developers?
– Sycorax
Commented Apr 15, 2014 at 5:18
• @Aksakal runif(1,-Inf, Inf) gives a warning and the documentation says that the limits must be finite. I prefer this to an error and this is consistend behavior in R. Commented Apr 15, 2014 at 10:54
• @user777, R developers are not real programmers. It's a waste of time to preach good coding practice to non programmers Commented Apr 15, 2014 at 11:43

Uniform distribution is not defined on unbounded domain. It has to be $[-\infty<a<x<b<\infty]$.

That was your main issue. The second issue is that when you integrate with MC you don't have to use uniform distribution. Your case is a good example of when not to use it, in fact. Look up importance sampling subject in any Monte Carlo tutorial. Obviously, normal or exponential distributions would be good places to start for your integrals.

In a nutshell, what you're doing is this: $\hat{I}=\sum_if(x_i)/p(x_i)$, where $p(x_i)$ is the PDF of your sampling distribution. In the most basic case it would uniform distribution, i.e. $\frac{1}{b-a}$. If you use normal distribution it'll be $\frac{e^{-x_i^2/2}}{\sqrt{2\pi}}$.

UPDATE: The MATLAB code would be like this:

x=randn(10000,1);
mean(exp(-x.^2)./normpdf(x))

ans =

1.7828

• So, I tried using rnorm and I get this: > f=function(n){ + x=rnorm(n) + g=exp(-x^2) + f=mean(g) + f + } > > f(10000) [1] 0.5800912 but using "integrate", I get this: > f <- function(x) {exp(-x^2)} > integrate(f,-Inf,Inf) 1.772454 with absolute error < 4.3e-06 What could I be doing wrong ? Commented Apr 14, 2014 at 16:01
• @integral, try using code markup so that we can read your code. Commented Apr 14, 2014 at 16:52
• Can't get code to work properly. Here is a pastebin link: pastebin.com/RKDLNbJa Commented Apr 14, 2014 at 17:57
• integral code markup works in questions and answers. Comments are not designed to hold chunks of code - if it's too long for the inline backquoted (inline-code) markup, put it at the bottom of your question via an edit. Commented Apr 15, 2014 at 1:31
• @integral, you didn't implement sampling right. Particularly, the change of measure part, i.e. where you divide by PDF. Commented Apr 15, 2014 at 2:12