Search code examples
rmontecarlo

Monte Carlo Integration of exp(-x^2) from -Infinity to +Infinity

I have to integrate ∫(-∞)^∞〖e^(-x^2 ) dx, I do not know how to tackle the limits and my current result is way off, this is my code: 

`ok<-function(x) exp(-x^2)
set.seed(999)
n.sim<-10^5
x1<-ok(runif(n.sim))
theta1<-mean(x1)
theta1
[1] 0.7468196 `
Solution

  • You can specify lower and upper limits for runif(). You need to make the limits large enough that it captures the tails of the distribution (the default (0,1) range is definitely not large enough!) but you need to keep it small enough for good coverage. +/- 3 seems reasonable (based on what we already know about the tails of the Normal distribution ...) You also have to multiply by the size of the domain ...

    set.seed(101)
mean(exp(-(runif(1e7,-3,3))^2))*6
## [1] 1.771848

    Compare with integrate():

    integrate(function(x)exp(-x^2), -3,3)
## 1.772415 with absolute error < 6.4e-11