4
$\begingroup$

Monte Carlo Estimates from Different distributions for $\theta = \int_{0}^1 \sin(x) dx$

So I understand that a monte carlo estimate for this integral should look like this $\hat\theta = \frac{1}{m} \sum_{i=1}^m sin(x_i)$, where $x_i \sim U(0,1)$

And in R I can see this works:

> m = 10^4
> x <- runif(m)
> theta.hat = mean(sin(x))
> cat(paste("Theta.hat = ", round(theta.hat, digits = 5)))
Theta.hat =  0.46435

Which is very close to the true value $1-\cos(1)\approx0.45969$

But I should also be able to estimate this by sampling from another distribution as well, say $Exp(1)$ where $f(x)=e^{-x}$

Since $\theta = \int_{0}^1 \frac{\sin(x)}{e^{-x}}\cdot e^{-x} dx= E(\frac{\sin(x)}{e^{-x}})$

Which I should have this estimator $\hat\theta = \frac{1}{m} \sum_{i=1}^m I(x_i<1) \frac{\sin(x_i)}{e^{-x_i}}$ with $x_i\sim Exp(1)$

When I try to do this in R however:

> m = 10^4
> accept = c()
> x <- rexp(m, rate = 1)
> for (i in x){
+   if(i<= 1){
+     accept=c(accept,(sin(i)/exp(-i)))
+   }
+ }
> mean(accept)
[1] 0.7255409

I get something that's significantly different from the actual value. I'm not sure what exactly it is that I'm doing wrong.

$\endgroup$

1 Answer 1

5
$\begingroup$

Your code in the latter case is taking the mean only over the cases where $x_i < 1$, which is incorrect. As in your equation, you need to take the mean over all the simulations. To do this, just replace mean(accept) with sum(accept)/m and you'll get a reasonable estimate.

A few other coding tips: You are presently working in Circle 2 of hell in the R inferno --- it is better to pre-construct your vectors to the desired length and then populate them rather than growing them element-by-element in this manner. Also, always set your seed for reproducibility when doing simulations involving PRNGs. Here is a better way to code the latter simulation:

#Generate estimator via importance sampling
m <- 10^4
set.seed(1)
x <- rexp(m, rate = 1)
SIM <- rep(0, m)
for (i in 1:m) { 
  if (x[i] < 1) {
    SIM[i] <- sin(x[i])/exp(-x[i]) } }

#Compare estimator to known value
THETA     <- 1-cos(1)
THETA.HAT <- mean(SIM)

THETA
[1] 0.4596977

THETA.HAT
[1] 0.4583018
$\endgroup$
1
  • $\begingroup$ Yes thanks. I actually just realized this myself when I decided to check the length of my accept variable. $\endgroup$ Mar 6, 2022 at 21:19

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.