Lower than expected coverage for importance sampling with simulation

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

π0f(x)dx=π01cos(x)2+x2dx

using the exponential distribution as the importance distribution

q(x)=λ expλx

and find the value of λ which gives the better approximation to the integral (it’s self-study). I recast the problem as the evaluation of the mean value μ of f(x) over [0,π]: the integral is then just πμ.

Thus, let p(x) be the pdf of XU(0,π), and let Yf(X): the goal now is to estimate

μ=E[Y]=E[f(X)]=Rf(x)p(x)dx=π01cos(x)2+x21π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 μ, and each time a checks is made on whether the 95% interval covers the actual mean or not.

As you can see, for λ=20 the actual coverage is just 0.19. And increasing B to values such as 106 doesn’t help (the coverage is even smaller, 0.15). Why is this happening?

Answer

Importance sampling is quite sensitive to the choice of importance distribution. Since you chose λ=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 π=3.14. So you want to use a λ that gives you such a range. I use λ=1.

enter image description here

Using λ=1 I will be able to explore the full integral space of 0 to π, and seems like only a few draws over π will be wasted. Now I rerun your code, and only change λ=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 λ, 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=104 to B=106, 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=104 is,
.19±1.96.19(1.19)100=.19±.0769=(.1131,.2669).

So you can’t really say that increasing B=106 significantly lowers the coverage probability.

In fact in your code for the same seed, change N=100 to N=1000, then with B=104, coverage probability is .123 and with B=106 coverage probability is .158.

Now, the confidence interval around .123 is
.123±1.96.123(1.123)1000=.123±.0203=(.102,.143).

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

Attribution
Source : Link , Question Author : DeltaIV , Answer Author : Greenparker

Leave a Comment