Example 4.6

We compare here the MC estimator for a tail probability estimation problem with the IS estimate.

import numpy as np
import matplotlib.pyplot as plt

xx = np.linspace(4, 20, 100000)

def p(x):
    return np.exp(-x**2/2)/np.sqrt(2*np.pi)

def q(x, mu, sigma):
    return np.exp(-(x-mu)**2/(2*sigma**2))/(np.sqrt(2*np.pi)*sigma)

def w(x, mu, sigma):
    return p(x)/q(x, mu, sigma)

I = np.trapz(p(xx), xx) # Numerical computation of the integral

print('Integral of p(x) from 4 to infinity: ', I)

N = 10000

x = np.random.normal(0, 1, N) # iid samples from p(x)

I_est_MC = (1/N) * np.sum(x > 4)
print('Monte Carlo estimate: ', I_est_MC)

mu = 6
sigma = 1

x_s = np.zeros(N)
weights = np.zeros(N)

for i in range(N):
    x_s[i] = np.random.normal(mu, sigma)
    weights[i] = w(x_s[i], mu, sigma)

I_est_IS = (1/N) * np.sum(weights * (x_s > 4))
print('Importance sampling estimate: ', I_est_IS)

Integral of p(x) from 4 to infinity:  3.16712429751607e-05
Monte Carlo estimate:  0.0001
Importance sampling estimate:  3.052061109150359e-05