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

- 算法步骤
输入：任意选定的建议分布（状态转移核）q，抽样的目标分布密度函数 p(x)，收敛步数 m，需要样本数 n 。
Step 1：随机选择一个初始值
，样本集合 =[] 。
Step 2：for t=0 to m+n：

文章插图
Step 3：返回样本集合。
一个简单的例子：

文章插图

# -*- coding:utf-8 -*-import randomimport numpy as npimport matplotlib.pyplot as pltdef mh(q, p, m, n):# randomize a numberx = random.uniform(0.1, 1)for t in range(0, m+n):x_sample = q.sample(x)try:accept_prob = min(1, p.prob(x_sample)*q.prob(x_sample, x)/(p.prob(x)*q.prob(x, x_sample)))except:accept_prob = 0u = random.uniform(0, 1)if u < accept_prob:x = x_sampleif t >= m:yield xclass Exponential(object):def __init__(self, scale):self.scale = scaleself.lam = 1.0 / scaledef prob(self, x):if x <= 0:raise Exception("The sample shouldn't be less than zero")result = self.lam * np.exp(-x * self.lam)return resultdef sample(self, num):sample = np.random.exponential(self.scale, num)return sample# 假设我们的目标概率密度函数p1(x)是指数概率密度函数scale = 5p1 = Exponential(scale)class Norm():def __init__(self, mean, std):self.mean = meanself.sigma = stddef prob(self, x):return np.exp(-(x - self.mean) ** 2 / (2 * self.sigma ** 2.0)) * 1.0 / (np.sqrt(2 * np.pi) * self.sigma)def sample(self, num):sample = np.random.normal(self.mean, self.sigma, size=num)return sample# 假设我们的目标概率密度函数p1(x)是均值方差分别为3,2的正态分布p2 = Norm(3, 2)class Transition():def __init__(self, sigma):self.sigma = sigmadef sample(self, cur_mean):cur_sample = np.random.normal(cur_mean, scale=self.sigma, size=1)[0]return cur_sampledef prob(self, mean, x):return np.exp(-(x-mean)**2/(2*self.sigma**2.0)) * 1.0/(np.sqrt(2 * np.pi)*self.sigma)# 假设我们的转移核方差为10的正态分布q = Transition(10)m = 100n = 100000 # 采样个数simulate_samples_p1 = [li for li in mh(q, p1, m, n)]plt.subplot(2,2,1)plt.hist(simulate_samples_p1, 100)plt.title("Simulated X ~ Exponential(1/5)")samples = p1.sample(n)plt.subplot(2,2,2)plt.hist(samples, 100)plt.title("True X ~ Exponential(1/5)")simulate_samples_p2 = [li for li in mh(q, p2, m, n)]plt.subplot(2,2,3)plt.hist(simulate_samples_p2, 50)plt.title("Simulated X ~ N(3,2)")samples = p2.sample(n)plt.subplot(2,2,4)plt.hist(samples, 50)plt.title("True X ~ N(3,2)")plt.suptitle("Transition Kernel N(0,10)simulation results")plt.show()

代码运行结果：

文章插图
多元情况
在很多情况下，我们的目标分布是多元联合概率分布
，我们可以用条件分布概率的比来计算联合概率的比，从而提升计算效率。即：

文章插图
其中：
并且
被称为满条件分布（full）。
而在用转移核对多元变量分布进行抽样的时候，可以对多元变量的每一个变量的条件分布一次分别进行抽样，从而实现对整个多元变量的一次抽样。
针对多元变量的情况，我们来讨论一下其中的一个特例，吉布斯采样（Gibbs ）。
Gibbs
吉布斯采样（Gibbs ）适用于多元变量联合分布的抽样和估计。比较特殊的地方是它的建议分布是当前变量
的满条件概率分布：
这是，接受概率 α=1：

文章插图
其中，
因为不难想象，假设对第 j 个变量采样前
，对其采完样后
，那么在 x′ 中分别去掉第 j 项，剩下的其余变量全部相同，所以它们的概率值也相同。