# 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. Apr 15 '14 at 1:40
• I already tried that, but I'm stuck at how to do a double integration in R, not on paper Apr 15 '14 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) Apr 15 '14 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.77192


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

(2*a) * sum(exp(-abs(runif(N, -a, a)))) / N
 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 ? Apr 15 '14 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)? Apr 15 '14 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
Apr 15 '14 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
Apr 15 '14 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. Apr 15 '14 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. Apr 14 '14 at 15:28
• @Aksakal Have you made this suggestion to the R developers?
– Sycorax
Apr 15 '14 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. Apr 15 '14 at 10:54
• @user777, R developers are not real programmers. It's a waste of time to preach good coding practice to non programmers Apr 15 '14 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)  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 ? Apr 14 '14 at 16:01
• @integral, try using code markup so that we can read your code. Apr 14 '14 at 16:52
• Can't get code to work properly. Here is a pastebin link: pastebin.com/RKDLNbJa Apr 14 '14 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. Apr 15 '14 at 1:31
• @integral, you didn't implement sampling right. Particularly, the change of measure part, i.e. where you divide by PDF. Apr 15 '14 at 2:12