3
$\begingroup$

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") 
}
$\endgroup$
  • $\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 Feb 6 '17 at 16:17
  • $\begingroup$ Ops, R code fixed! $\endgroup$ – Paul Feb 6 '17 at 16:24
  • $\begingroup$ You might find this helpful: stats.stackexchange.com/questions/250934/… $\endgroup$ – Taylor Feb 6 '17 at 16:46
  • $\begingroup$ @Taylor Nope, I can't get it :-) I edited my initial post. $\endgroup$ – Paul Feb 6 '17 at 22:02
1
$\begingroup$

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
$\endgroup$
  • $\begingroup$ Nooo you beat me to it :) your me code is equivalent to what I did here right? $\endgroup$ – DeltaIV Feb 11 '17 at 9:25
  • 1
    $\begingroup$ Yes it is the same notion. $\endgroup$ – Xi'an Feb 11 '17 at 9:28
0
$\begingroup$

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.

$\endgroup$
  • $\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 Feb 10 '17 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 Feb 10 '17 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 Feb 10 '17 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 Feb 10 '17 at 17:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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