An unnormalized distribution assigns density weights up to an unknown constant. Writing p(x) = p̃(x)/Zp and q(x) = q̃(x)/Zq separates each density into its unnormalized form and normalization constant. section 2.7
Importance sampling uses this representation to avoid requiring known partition functions:
- Compute unnormalized importance weights from ratios such as p̃(x)/q̃(x). section 2.7
- Estimate ratios of normalization constants through weighted samples. section 2.7
- Analyze estimator behavior using the weak law of large numbers, variance, and Chebyshev’s inequality. section 2.7
Importance-sampling implementation
A target density can be implemented as an unnormalized function such as x^(k−1) exp(−x²/2), while samples come from a scaled Gaussian proposal. The estimator divides target values by proposal values and uses the resulting weights. 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 complete importance-sampling implementation, including target and proposal densities, sampling, and weight computation. section 2.7
