Metropolis-Hastings sampling constructs a Markov chain whose stationary distribution matches a target distribution p(x), such as a posterior or density. From the current state x, it proposes x′ from q(x′|x), then accepts or rejects the proposal using an MH ratio that preserves detailed balance. section 2.6 section 2.9
The method proceeds by:
- Selecting a target distribution p(x) and a proposal distribution q(x′|x). section 2.6
- Proposing a new state x′ conditional on the current state x. section 2.6
- Computing the MH ratio, derived as the correction required for detailed balance. section 2.6
- Accepting x′ with probability r(x′|x), or retaining the current state with probability 1 − r(x′|x). section 2.6
Sampling a Gaussian mixture
A from-scratch implementation defines target_pdf, proposal_pdf, and sample methods in an mh_gauss class. It samples a multivariate Gaussian-mixture target using a multivariate normal proposal, producing samples resembling two two-dimensional Gaussian components with means [4,0] and [−4,0], covariances [[2,1],[1,1]] and [[1,0],[0,1]], and mixture proportions [0.4,0.6]. section 2.6
Figure 2.13 Metropolis-Hastings pseudo-code
Use the pseudocode to follow the proposal, ratio, acceptance, and rejection loop. section 2.6
Listing 2.5 Metropolis-Hastings sampling
import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import uniform
from scipy.stats import multivariate_normal
np.random.seed(42)
class mh_gauss:
def __init__(self, dim, K, num_samples, target_mu, target_sigma,
target_pi, proposal_mu, proposal_sigma):
self.dim = dim
self.K = K #1
self.num_samples = num_samples #1
self.target_mu = target_mu #1
self.target_sigma = target_sigma #1
self.target_pi = target_pi #1
self.proposal_mu = proposal_mu #1
self.proposal_sigma = proposal_sigma #2
self.n_accept = 0
self.alpha = np.zeros(self.num_samples) #3
self.mh_samples = np.zeros((self.num_samples, #3
self.dim)) #3
def target_pdf(self, x):
prob = 0
for k in range(self.K): #4
prob += self.target_pi[k]*\ #4
multivariate_normal.pdf(x,self #4
.target_mu[:,k],self.target_sigma[:,:,k]) #4
#end for #4
return prob #4
def proposal_pdf(self, x):
return multivariate_normal.pdf(x, self
.proposal_mu, self.proposal_sigma) #5
def sample(self):
x_init = multivariate_normal
.rvs(self.proposal_mu, self.proposal_sigma, 1) #6
self.mh_samples[0,:] = x_init
for i in range(self.num_samples-1):
x_curr = self.mh_samples[i,:]
x_new = multivariate_normal.rvs(x_curr, self
.proposal_sigma, 1)
self.alpha[i] = self.proposal_pdf(x_new) /
self.proposal_pdf(x_curr) #7
self.alpha[i] = self.alpha[i] * #7
(self.target_pdf(x_new)/self.target_pdf(x_curr)) #7
r = min(1, self.alpha[i]) #8
u = uniform.rvs(loc=0, scale=1, size=1)
if (u <= r):
self.n_accept += 1
self.mh_samples[i+1,:] = x_new #9
else:
self.mh_samples[i+1,:] = x_curr #10
#end for
print("MH acceptance ratio: ", self.n_accept/float(self.num_samples))
if __name__ == "__main__":
dim = 2
K = 2
num_samples = 5000
target_mu = np.zeros((dim, K))
target_mu[:,0] = [4,0]
target_mu[:,1] = [-4,0]
target_sigma = np.zeros((dim, dim, K))
target_sigma[:,:,0] = [[2,1],[1,1]]
target_sigma[:,:,1] = [[1,0],[0,1]]
target_pi = np.array([0.4, 0.6])
proposal_mu = np.zeros((dim,1)).flatten()
proposal_sigma = 10*np.eye(dim)
mhg = mh_gauss(dim, K, num_samples, target_mu, target_sigma, target_pi,
proposal_mu, proposal_sigma)
mhg.sample()
plt.figure()
plt.scatter(mhg.mh_samples[:,0], mhg.mh_samples[:,1], label='MH samples')
plt.grid(True); plt.legend()
plt.title("Metropolis-Hastings Sampling of 2D Gaussian Mixture")
plt.xlabel("X1"); plt.ylabel("X2")
plt.show()
Use the implementation listing to inspect target evaluation, proposal evaluation, sampling, and acceptance tracking. section 2.6
