Evaluate integral with Importance sampling method in R

Question

I need to evaluate this integral with importance sampling method $$\int_{0}^{\pi}\frac{1}{\cos(x)^2+x^2}dx$$ using this weight function (exponential) $$f(x)=\lambda\exp\{-\lambda x\}$$

This is my attempt

set.seed(1)
f <- function(x){
    1 / (cos(x)^2+x^2)
}
mc <- function(lambda, f, B){
    x <- rexp(B, lambda)
    f(x[x<pi]) / dexp(x[x<pi], lambda)
}
lambda <- 1
B <- 10000
I <- mc(lambda, f, B)
mean(I)
1.648396

How can I find which $\lambda$ minimize approximation error?

First edit

I try with this but the value that give me the lower variance is not giving the right mean.

lambda <- seq(0.1, 3, 0.1)
for(i in lambda){
    I <- mc(i, f, B)
    cat("lambda = ", i, "mean = ", mean(I), "var = ", var(I), "\n") 
}

Answer 1

I think there is a fundamental mistake in your code:

mc <- function(lambda, f, B){
    x <- rexp(B, lambda)
    f(x[x<pi]) / dexp(x[x<pi], lambda)
}

since it produces a sample of an exponential truncated to $(0,\pi)$, with a random size. Therefore, the importance weight should be the inverse of the density of the truncated Exponential,$$\lambda\exp\{-\lambda x\}\big/1-\exp\{-\lambda\pi\}$$that is

md <- function(lambda, f, B){
    x <- rexp(B, lambda)
    pexp(pi,lambda)*f(x[x<pi]) / dexp(x[x<pi], lambda)
}

Comparing the outputs shows why one md works and the other mc does not:

> integrate(f,0,pi)$val
[1] 1.581188
> lambda=.1;mean(mc(lambda,f,B));mean(md(lambda,f,B))
[1] 5.800769
[1] 1.556356
> lambda=.5;mean(mc(lambda,f,B));mean(md(lambda,f,B))
[1] 1.991859
[1] 1.588999
> lambda=2;mean(mc(lambda,f,B));mean(md(lambda,f,B))
[1] 1.566827
[1] 1.597629

The alternative to the truncation is to return the entire Exponential sample while taking $f(x)=0$ for $x>\pi$:

me <- function(lambda, f, B){
    x <- rexp(B, lambda)
    f(x)*(x<pi) / dexp(x, lambda)
}

> lambda=.1;mean(mc(lambda,f,B));mean(md(lambda,f,B));mean(me(lambda,f,B));
[1] 5.869764
[1] 1.575606
[1] 1.574132

This being fixed, looking for the optimal value of $\lambda$ can be operated by minimising the Monte Carlo variance over $\lambda$ with the very same Exponential sample, in order to avoid Monte Carlo variability

expone=rexp(B,1) #standard Exponential
target <- function(lambda,expone,B){
  explambda=expone/lambda
  fexplambda=f(explambda)*(explambda<pi)/dexp(explambda,lambda)
  return(var(fexplambda))}

which leads to

> optimise(target,c(.01,5),expone=expone,B=1e6)
$minimum
[1] 0.9462493

$objective
[1] 0.2354651

Answer 2

Your mean and standard deviation are random. You will never have exactly the right mean (with probability one). Notice how if you run the last portion again, all the numbers change. Maybe you would like a confidence interval: "mean" $\pm$ $2\sqrt{\text{"var"}}$.

Also, disregard my comment, because that applies to self-normalized importance sampling.