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

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)

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") 
  • $\begingroup$ B is not declared in your code, so the code doesn't run. Please fix your code (it's presumably the number of MC samples). $\endgroup$
    – DeltaIV
    Commented Feb 6, 2017 at 16:17
  • $\begingroup$ Ops, R code fixed! $\endgroup$
    – Paul
    Commented Feb 6, 2017 at 16:24
  • $\begingroup$ You might find this helpful: stats.stackexchange.com/questions/250934/… $\endgroup$
    – Taylor
    Commented Feb 6, 2017 at 16:46
  • $\begingroup$ @Taylor Nope, I can't get it :-) I edited my initial post. $\endgroup$
    – Paul
    Commented Feb 6, 2017 at 22:02

2 Answers 2


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){

which leads to

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

[1] 0.2354651
  • $\begingroup$ Nooo you beat me to it :) your me code is equivalent to what I did here right? $\endgroup$
    – DeltaIV
    Commented Feb 11, 2017 at 9:25
  • 1
    $\begingroup$ Yes it is the same notion. $\endgroup$
    – Xi'an
    Commented Feb 11, 2017 at 9:28

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.

  • $\begingroup$ I don't think this is the right way to go, but even if it was, the confidence interval for the sample mean would be $\mu\pm 2 \frac{\sigma}{\sqrt{B}}$. You have to divide by the (square root of) sample size. $\endgroup$
    – DeltaIV
    Commented Feb 10, 2017 at 10:03
  • $\begingroup$ @DeltaIV $\frac{\sigma}{\sqrt{B}} = \sqrt{\text{Var}(\bar{X})}$. I am applying the law of large numbers/CLT. When Paul says he isn't getting the right mean, I assume he's wondering why his answers don't correspond exactly with his known answer. $\endgroup$
    – Taylor
    Commented Feb 10, 2017 at 16:41
  • $\begingroup$ ok, you didn't specify what variance you were talking about, so I thought you were talking about the sample variance of the Importance Sampling sample. I think Paul is actually wondering about why the $\lambda$ which corresponds to minimal sample variance in its code (minimal sample variance is what you usually aim for, when selecting the importance distribution) doesn't also minimize the error (difference with the exact value of the integral). I'll elaborate in an answer in the next days. $\endgroup$
    – DeltaIV
    Commented Feb 10, 2017 at 16:50
  • $\begingroup$ Right, and I'm saying the difference between the exact expectation (constant) and the sample mean (random variable) is random. So it shouldn't necessarily be the smallest difference among all his attempts. The way to go would be to ask, "how should I space my $\lambda$ grid in order to get the best possible estimate for the optimal instrumental density." That would be a cool question $\endgroup$
    – Taylor
    Commented Feb 10, 2017 at 17:10

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.