Importance sampling estimates an integral or expectation under a target distribution p(x) by drawing samples from an easier proposal distribution q(x). It concentrates samples in regions where both p(x) and |f(x)| are large, then corrects the proposal with the importance weight w(x)=p(x)/q(x). section 2.7
The basic estimator follows this sequence:
- Choose a proposal q(x) that is easier to sample and represents the relevant regions of the target and integrand. section 2.7
- Draw N samples xᵢ from q(x), compute wᵢ=p(xᵢ)/q(xᵢ), and weight each f(xᵢ) by wᵢ. section 2.7
- Approximate the Monte Carlo integral by averaging the weighted values: (1/N)Σᵢ wᵢf(xᵢ). section 2.7
Noninteger-degree-of-freedom chi target
For f(x)=2 sin(πx/1.5) with x≥0, a proposal based on a Gaussian is sampled, restricted to positive values, and compared with an unnormalized chi target. The proposal density is scaled by c so that p(x)<c q(x) for all x≥0. The estimate is compared with a ground truth obtained by numerical integration as the sample count increases. section 2.7
Figure 2.15 Importance sampler pseudo-code
Use the pseudo-code to trace proposal sampling, positivity filtering, weight computation, and weighted averaging. section 2.7
Listing 2.6 Importance sampling
import numpy as np
import matplotlib.pyplot as plt
from scipy.integrate import quad
from scipy.stats import multivariate_normal
np.random.seed(42)
class importance_sampler:
def __init__(self, k=1.5, mu=0.8, sigma=np.sqrt(1.5), c=3):
self.k = k #1
self.mu = mu #2
self.sigma = sigma #2
self.c = c #3
def target_pdf(self, x):
return (x**(self.k-1)) * np.exp(-x**2/2.0) #4
def proposal_pdf(self, x):
return self.c * 1.0/np.sqrt(2*np.pi*1.5) *
np.exp(-(x-self.mu)**2/(2*self.sigma**2)) #5
def fx(self, x):
return 2*np.sin((np.pi/1.5)*x) #6
def sample(self, num_samples):
x = multivariate_normal.rvs(self.mu,
self.sigma, num_samples) #7
idx = np.where(x >= 0)
x_pos = x[idx] #8
isw = self.target_pdf(x_pos) /
self.proposal_pdf(x_pos) #9
fw = (isw/np.sum(isw))*self.fx(x_pos)
f_est = np.sum(fw) #10
return isw, f_est
if __name__ == "__main__":
num_samples = [10, 100, 1000, 10000, 100000, 1000000]
F_est_iter, IS_weights_var_iter = [], []
for k in num_samples:
IS = importance_sampler()
IS_weights, F_est = IS.sample(k)
IS_weights_var = np.var(IS_weights/np.sum(IS_weights))
F_est_iter.append(F_est)
IS_weights_var_iter.append(IS_weights_var)
#ground truth (numerical integration)
k = 1.5
I_gt, _ = quad(lambda x: 2.0*np.sin((np.pi/1.5)*x)*(x**(k-1))*np
.exp(-x**2/2.0), 0, 5)
#generate plots
plt.figure()
xx = np.linspace(0,8,100)
plt.plot(xx, IS.target_pdf(xx), '-r', label='target pdf p(x)')
plt.plot(xx, IS.proposal_pdf(xx), '-b', label='proposal pdf q(x)')
plt.plot(xx, IS.fx(xx) * IS.target_pdf(xx), '-k', label='p(x)f(x) integrand')
plt.grid(True); plt.legend(); plt.xlabel("X1"); plt.ylabel("X2")
plt.title("Importance Sampling Components")
plt.show()
plt.figure()
plt.hist(IS_weights, label = "IS weights")
plt.grid(True); plt.legend();
plt.title("Importance Weights Histogram")
plt.show()
plt.figure()
plt.semilogx(num_samples, F_est_iter, label = "IS Estimate of E[f(x)]")
plt.semilogx(num_samples, I_gt*np.ones(len(num_samples)), label = "Ground
Truth")
plt.grid(True); plt.legend(); plt.xlabel('iterations'); plt.ylabel("E[f(x)]
estimate")
plt.title("IS Estimate of E[f(x)]")
plt.show()
Inspect the from-scratch Python implementation defining target_pdf, proposal_pdf, fx, and sample. section 2.7
