Updated 2020-09-23

Source

Importance Sampling : Example 1

Suppose you want to evaluate the integral:

\[ \int_a^b h(x) dx = \int_{0}^{10} \exp(-2 | x- 5| ) dx \] The simplest way is to generate \(U_i \sim \text{Unif}(0,10)\) and look at the sample mean \((b-a) \times h(U_i)\).

Importance Sampling

\[ \begin{align*} \int_a^b h(x) dx & = \int_{0}^{10} \exp(-2 | x- 5| ) dx \\ & = \int_{0}^{10} 10 \exp(-2 |x- 5| ) \frac{1}{10} dx \\ & = E (10\exp(-2 \lvert U - 5 \rvert )); \; U \sim \text{Unif}(0,10) \end{align*} \]

R codes

set.seed(123)
U <-  runif(1e4, 0, 10)
Y <- 10* exp(-2*abs(U - 5))
(I <- mean(Y))
## [1] 1.01854
h <- function(u) exp(-2*abs(u - 5))
integrate(h, 0,10)$val
## [1] 0.9999546
  • They are close !

Error

  • How do you estimate the error ?
  • Monte Carlo Error: use the variance of generated samples.
(Error <- var(Y))
## [1] 4.040725
  • Now this is quite large !

Improving the accuracy

  1. This method has a high variance, or, estimation error.
  2. The function h in this case is peaked at 5, and decays quickly elsewhere.
  3. Under the \(U[0,10]\) density, many of the points are contributing very little to this expectation.
  4. Something more like a gaussian function with a peak at 5 and small variance, (e.g. \(\sigma^2 = 1\)) would provide greater precision.

Better

Importance Sampling

Write the same integral as an expectation of a function under Normal distribution:

\[ \begin{align*} I & = \int_{0}^{10} 10 \exp(-2 | x- 5| ) \frac{\frac{1}{10}}{\phi(x-5)} \phi(x-5) dx \\ & = \int_{\mathbb{R}} 10 \exp(-2 | x- 5| ) \frac{\frac{1}{10} I(0 \le x \le 10)}{\phi(x-5)} \phi(x-5) dx \\ & \text{where } \phi(x-5) = \frac{1}{\sqrt{2\pi}} e^{-(x-5)^2/2} \\ & = E \{ h(X)w(X) \}, \text{ for } X \sim N(5,1) \end{align*} \] Q. What are \(h(x)\) and \(w(x)\) in the above example?

Importance Sampling

\[ \begin{align*} I & = \int_{\mathbb{R}} 10 \exp(-2 | x- 5| ) \frac{\frac{1}{10} I(0 \le x \le 10)}{\phi(x-5)} \phi(x-5) dx \\ & = E \{ h(X)w(X) \}, \text{ for } X \sim N(5,1) \end{align*} \]

  • We can compare sides and find \(h(x)\) and \(w(x)\):

\[ \begin{align*} & \Rightarrow h(x) = 10 \exp(-2 | x- 5| ) \\ & \Rightarrow w(x) = \frac{\text{density function of } U(0,10)}{\text{density function of } N(5,1)} \\ & \text{in R-language, }w(x) = \frac{\text{dunif(x, 0, 10)}}{\text{dnorm(x, 5, 1)}} \end{align*} \]

R code

w <- function(x) dunif(x, 0, 10)/dnorm(x, mean=5, sd=1)
f <- function(x) 10*exp(-2*abs(x-5))
X=rnorm(1e5,mean=5,sd=1)
Y=w(X)*f(X)

(I <- mean(Y))
## [1] 0.9983512
(Error <- var(Y))
## [1] 0.3552715

Same integral but the new error is 1/10th of the previous error.

Example 2

The standard double exponential or Laplace density is given by:

\[ p(x) = \frac{1}{2} e^{-\lvert x \rvert} \]

Calculate the second moment \(E(X^2)\) for Laplace using the Normal(0,4) density as the proposal. In other words, calculate the integral:

\[ \int_{-\infty}^{\infty} \frac{1}{2} x^2 e^{-\lvert x \rvert} dx \]

Report both your estimate and the variance of your estimate.

Could you also calculate the true value?

Solution

X <- rnorm(1e5, sd=4)
Y <- (X^2) * .5 * exp(-abs(X))/dnorm(X, sd=4)
(I <- mean(Y))
## [1] 1.997235
(Error <- var(Y))
## [1] 0.8958159

Think: What happens if we use \(N(0,1)\) instead? Or a Cauchy instead?

True Value

The true value of the integral is 2. You can find that by integrating by hand, or, evaluating it in R.

h <- function(x) 0.5*x^2*exp(-abs(x))
integrate(h,-Inf,Inf)$val
## [1] 2