如何快速理解马尔科夫链蒙特卡洛法？( 五 ) _概率分布

同理：
我们可以发现因为 α 接受分布永远等于 1，那么吉布斯采样是不会拒绝样本的。
这里需要注意的是，这里的目标分布的单个元的分布至少应该是可以采样的，当然如果单个元很复杂的话，还是可以再用一次 M-H 的方法来为这个单元的分布进行采样（存疑）。
吉布斯采样算法
输入：抽样的目标分布密度函数 p(x)，转移核
收敛步数 m，需要样本数 n 。
Step 1：随机选择一个样本
，样本集合 =[]
Step 2：
for i=1 to m+n:
for j=1 to k:
根据：
抽取
，并赋值给
。
if (i >= m) 将
加入到中。
Step 3：返回。
来看个例子吧。假设我们有二维正态分布：
【如何快速理解马尔科夫链蒙特卡洛法？】

文章插图
假设我们期望抽样的二元正态分布是 μ=(5,8)，协方差矩阵为：

文章插图

import matplotlib.pyplot as pltfrom mpl_toolkits.mplot3d import Axes3Dimport numpy as npclass Transition():def __init__(self, mean, cov):self.mean = meanself.sigmas = []for i in range(K):self.sigmas.append(np.sqrt(cov[i][i]))self.rho = cov[0][1]/(self.sigmas[0] * self.sigmas[1])def sample(self, id1, id2_list, x2_list):id2 = id2_list[0]# only consider two dimensionx2 = x2_list[0]# only consider two dimensioncur_mean = self.mean[id1] + self.rho*self.sigmas[id1]/self.sigmas[id2] * (x2-self.mean[id2])cur_sigma = (1-self.rho**2) * self.sigmas[id1]**2return np.random.normal(cur_mean, scale=cur_sigma, size=1)[0]def gibbs(p, m, n):# randomize a numberx = np.random.rand(K)for t in range(0, m+n):for j in range(K):total_indexes = list(range(K))total_indexes.remove(j)left_x = x[total_indexes]x[j] = p.sample(j, total_indexes, left_x)if t >= m:yield xmean = [5, 8]cov = [[1, 0.5], [0.5, 1]]K = len(mean)q = Transition(mean, cov)m = 100n = 1000gib = gibbs(q, m, n)simulated_samples = []x_samples = []y_samples = []for li in gib:x_samples.append(li[0])y_samples.append(li[1])fig = plt.figure()ax = fig.add_subplot(131, projection='3d')hist, xedges, yedges = np.histogram2d(x_samples, y_samples, bins=100, range=[[0,10],[0,16]])xpos, ypos = np.meshgrid(xedges[:-1], yedges[:-1])xpos = xpos.ravel()ypos = ypos.ravel()zpos = 0dx = xedges[1] - xedges[0]dy = yedges[1] - yedges[0]dz = hist.flatten()ax.bar3d(xpos, ypos, zpos, dx, dy, dz, zsort='average')ax = fig.add_subplot(132)ax.hist(x_samples, bins=50)ax.set_title("Simulated on dim1")ax = fig.add_subplot(133)ax.hist(y_samples, bins=50)ax.set_title("Simulated on dim2")plt.show()

代码运行结果：

文章插图
从马尔科夫链蒙特卡洛法中我们得到的样本序列，相邻点是相关的，因此如果需要独立样本，我们可以在该样本序列中再次进行随机抽样，比如每隔一段时间抽取一次样本，将这样得到的样本集合作为独立样本集合。
另外，关于“燃烧期”，即过多少步后马尔科夫链是趋向稳定分布的，通常是经验性的。比如书上说，可以游走一段时间后采集的样本均值和过一段时间后采集的样本均值做个比较，如果均值稳定，那可以认为马尔科夫链收敛。
总结
本文回顾了马尔科夫链蒙特卡洛法的相关基本知识点、-, Gibbs基本原理。下一篇文章，我们一起来看看 MCMC 大法在实际应用中的一些例子。
[1] 李航. 统计学习方法第二版[J]. 2019.
[2]