9
$\begingroup$

I was trying to answer the question Evaluate integral with Importance sampling method in R. Basically, the user needs to calculate

$$\int_{0}^{\pi}f(x)dx=\int_{0}^{\pi}\frac{1}{\cos(x)^2+x^2}dx$$

using the exponential distribution as the importance distribution

$$q(x)=\lambda\ \exp^{-\lambda x}$$

and find the value of $\lambda$ which gives the better approximation to the integral (it's self-study). I recast the problem as the evaluation of the mean value $\mu$ of $f(x)$ over $[0,\pi]$: the integral is then just $\pi\mu$.

Thus, let $p(x)$ be the pdf of $X\sim\mathcal{U}(0,\pi)$, and let $Y\sim f(X)$: the goal now is to estimate

$$\mu=\mathbb{E}[Y]=\mathbb{E}[f(X)]=\int_{\mathbb{R}}f(x)p(x)dx=\int_{0}^{\pi}\frac{1}{\cos(x)^2+x^2}\frac{1}{\pi}dx$$

using importance sampling. I performed a simulation in R:

# clear the environment and set the seed for reproducibility
rm(list=ls())
gc()
graphics.off()
set.seed(1)

# function to be integrated
f <- function(x){
    1 / (cos(x)^2+x^2)
}

# importance sampling
importance.sampling <- function(lambda, f, B){
    x <- rexp(B, lambda) 
    f(x) / dexp(x, lambda)*dunif(x, 0, pi)
}

# mean value of f
mu.num <- integrate(f,0,pi)$value/pi

# initialize code
means  <- 0
sigmas <- 0
error  <- 0
CI.min <- 0
CI.max <- 0
CI.covers.parameter <- FALSE

# set a value for lambda: we will repeat importance sampling N times to verify
# coverage
N <- 100
lambda <- rep(20,N)

# set the sample size for importance sampling
B <- 10^4

# - estimate the mean value of f using importance sampling, N times
# - compute a confidence interval for the mean each time
# - CI.covers.parameter is set to TRUE if the estimated confidence 
#   interval contains the mean value computed by integrate, otherwise
# is set to FALSE
j <- 0
for(i in lambda){
    I <- importance.sampling(i, f, B)
    j <- j + 1
    mu <- mean(I)
    std <- sd(I)
    lower.CB <- mu - 1.96*std/sqrt(B)  
    upper.CB <- mu + 1.96*std/sqrt(B)  
    means[j] <- mu
    sigmas[j] <- std
    error[j] <- abs(mu-mu.num)
    CI.min[j] <- lower.CB
    CI.max[j] <- upper.CB
    CI.covers.parameter[j] <- lower.CB < mu.num & mu.num < upper.CB
}

# build a dataframe in case you want to have a look at the results for each run
df <- data.frame(lambda, means, sigmas, error, CI.min, CI.max, CI.covers.parameter)

# so, what's the coverage?
mean(CI.covers.parameter)
# [1] 0.19

The code is basically a straightforward implementation of importance sampling, following the notation used here. The importance sampling is then repeated $N$ times to get multiple estimates of $\mu$, and each time a checks is made on whether the 95% interval covers the actual mean or not.

As you can see, for $\lambda=20$ the actual coverage is just 0.19. And increasing $B$ to values such as $10^6$ doesn't help (the coverage is even smaller, 0.15). Why is this happening?

$\endgroup$
2
  • 1
    $\begingroup$ Using an infinite support importance function for a finite support integral is not optimal as a part of the simulations is used to simulate zeros, so to speak. At least truncate the exponential at $\pi$, which is is easy to do and simulate. $\endgroup$
    – Xi'an
    Feb 11, 2017 at 8:54
  • $\begingroup$ @Xi'an sure, I agree, if I had to evaluate that integral by Importance Sampling, I wouldn't use that importance distribution, but I was trying to answer the original question, which required using the exponential distribution. My problem was that, even if this approach is far from optimal, coverage should still increase (on average) as $B\to\infty$. And that's what Greenparker showed. $\endgroup$
    – DeltaIV
    Feb 11, 2017 at 14:40

1 Answer 1

3
+100
$\begingroup$

Importance sampling is quite sensitive to the choice of importance distribution. Since you chose $\lambda = 20$, the samples that you draw using rexp will have a mean of $1/20$ with variance $1/400$. This is the distribution you get

enter image description here

However, the integral you want to evaluate goes from 0 to $\pi =3.14$. So you want to use a $\lambda$ that gives you such a range. I use $\lambda = 1$.

enter image description here

Using $\lambda = 1$ I will be able to explore the full integral space of 0 to $\pi$, and seems like only a few draws over $\pi$ will be wasted. Now I rerun your code, and only change $\lambda = 1$.

# clear the environment and set the seed for reproducibility
rm(list=ls())
gc()
graphics.off()
set.seed(1)

# function to be integrated
f <- function(x){
  1 / (cos(x)^2+x^2)
}

# importance sampling
importance.sampling <- function(lambda, f, B){
  x <- rexp(B, lambda) 
  f(x) / dexp(x, lambda)*dunif(x, 0, pi)
}

# mean value of f
mu.num <- integrate(f,0,pi)$value/pi

# initialize code
means  <- 0
sigmas <- 0
error  <- 0
CI.min <- 0
CI.max <- 0
CI.covers.parameter <- FALSE

# set a value for lambda: we will repeat importance sampling N times to verify
# coverage
N <- 100
lambda <- rep(1,N)

# set the sample size for importance sampling
B <- 10^4

# - estimate the mean value of f using importance sampling, N times
# - compute a confidence interval for the mean each time
# - CI.covers.parameter is set to TRUE if the estimated confidence 
#   interval contains the mean value computed by integrate, otherwise
# is set to FALSE
j <- 0
for(i in lambda){
  I <- importance.sampling(i, f, B)
  j <- j + 1
  mu <- mean(I)
  std <- sd(I)
  lower.CB <- mu - 1.96*std/sqrt(B)  
  upper.CB <- mu + 1.96*std/sqrt(B)  
  means[j] <- mu
  sigmas[j] <- std
  error[j] <- abs(mu-mu.num)
  CI.min[j] <- lower.CB
  CI.max[j] <- upper.CB
  CI.covers.parameter[j] <- lower.CB < mu.num & mu.num < upper.CB
}

# build a dataframe in case you want to have a look at the results for each run
df <- data.frame(lambda, means, sigmas, error, CI.min, CI.max, CI.covers.parameter)

# so, what's the coverage?
mean(CI.covers.parameter)
#[1] .95

If you play around with $\lambda$, you will see that if you make it really small (.00001) or large, the coverage probabilities will be bad.

EDIT-------

Regarding the coverage probability decreasing once you go from $B = 10^4$ to $B = 10^6$, that is just a random occurrence, based on the fact that you use $N = 100$ replications. The confidence interval for the coverage probability at $B = 10^4$ is, $$.19 \pm 1.96*\sqrt{\dfrac{.19*(1-.19)}{100}} = .19 \pm .0769 = (.1131, .2669)\,.$$

So you can't really say that increasing $B = 10^6$ significantly lowers the coverage probability.

In fact in your code for the same seed, change $N = 100$ to $N = 1000$, then with $B = 10^4$, coverage probability is .123 and with $B = 10^6$ coverage probability is $.158$.

Now, the confidence interval around .123 is $$.123 \pm 1.96\sqrt{\dfrac{.123*(1 - .123)}{1000}} = .123 \pm .0203 = (.102, .143)\,. $$

Thus, now with $N = 1000$ replications, you get that the coverage probabiulity significantly increases.

$\endgroup$
3
  • $\begingroup$ Yes, I know that the coverage change with $\lambda$: in particular the best coverage is obtained for $0.1<\lambda<2$. Now, I understand that since the CI for the sample mean is based on CLT, it's an asymptotic result. Thus, it may well be that changing $\lambda$ influences the number of samples needed to approach the "asymptotic regime", so to speak. But the point is, why with $\lambda =20$ the coverage decreases from sample size $10^4$ to sample size $10^6$? Surely it should increase, if poor coverage was only due to an high $\lambda$ value? $\endgroup$
    – DeltaIV
    Feb 10, 2017 at 16:41
  • 1
    $\begingroup$ @DeltaIV I made an edit to answer this question. The gist is, $N = 100$ is not enough replications to say anything with certainty. $\endgroup$ Feb 10, 2017 at 16:55
  • 1
    $\begingroup$ ah brilliant! I didn't think of forming the confidence interval for the coverage proportion itself, instead than just for the mean. Just as a nitpick, I wouldn't have used the Wald confidence interval for the confidence interval of a proportion. However, since the proportion is away from 0 and 1 and the number of replicates is (in your second case, $N=1000$) relatively large, probably using the Wilson or Jeffreys interval wouldn't have made any difference. I will wait just a bit to see if there are other answers, but I'd say you fully deserve the +100 :) $\endgroup$
    – DeltaIV
    Feb 10, 2017 at 17:22

Your Answer

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

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