分布是一系列数字的规律组合。如果在收集了历史中的几百个数据后，我想知道这群数据背后的发射机制是什么，那么就得去寻找这个分布。当然这里的重点不是寻找分布，而是在已知分布的情况下，如何模拟这个机制发射出来的一系列数字呢？

MCMC（Markov Chain Monte Carlo）是马尔科夫链下的蒙特卡洛方法，因为马尔科夫链在满足某些条件下具有平稳分布，如果能够将平稳分布与目标分布联系起来，那么就可以达到对目标分布进行抽样的目的。这里主要介绍的是Metropolis Hasting 算法和Gibbs sampling 算法。

一、Metropolis Hasting

1、算法理解

我们的目标是对Target Distribution进行抽样，首先，我们要引入一条具有平稳分布的马氏链，这条马氏链收敛的平稳分布我们称为Proposal Distribution，而这条马氏链的表现形式是概率转移矩阵

，状态空间

，状态空间也即是Proposal distribution的所有可能取值集合。

如何根据这条马氏链求得目标分布呢？这里由马氏链的细致平稳性引入。

是目标分布下的随机变量，

是proposal distribution下的随机变量。

（1.1）成立。（由马氏链的细致平稳性得到，表示i，j状态之间的能量转换相等）

（1.2）（因为

与

是两个不同的分布）

（1.3）

为了使（1.2）式成立，所以引入了接受率

。其中

即将不等式的左右两边互相相乘，即可得到式子（1.3）。接受率

表示是否决定抽取下一个样本（i.e., 接受样本j），因此我们需要将这个概率实现，因为在实际抽样过程中，决定抽样和不抽样是一个二元过程，而不是说以多大的概率决定抽样。这个概率实现可以用伯努利分布，也可以用均匀分布：当均匀分布下的数值小于接受率时，决定抽样，反之不抽样。

以上就是Metropolis抽样方法的全部内容了，而Metropolis hasting 算法则对接受率做了一点改进。当接受率太小的时候，我们很难从当前的样本值跳到其他状态，所以对

进行了扩大。将

中的较大值扩充到1（即一定会抽取下个样本），另外一个值等比例扩大。经过计算可以得到表达式

。在计算接受率的过程中，我们就会发现，目标分布的常数项被抵消了，也去除了归一化的过程。

2、proposal distribution的选取

当proposal distribution与目标分布越靠近时，抽取的样本也就越合理。但是proposal distribution下的马氏链如何确定，两个分布的距离如何衡量，这些也都是可以继续探讨也需要权衡的问题。

3、共轭的正态分布示例

已知，

未知，在贝叶斯统计下，

是一个随机变量，其先验分布为

已知。如何利用Metropolis-Hasting算法，在观察数据Y下求得

后验分布得期望和方差？

我们用M-H抽样算法来检验上面得后验分布是否准确。即在已知得各参数和观测值y下抽出一系列的

。

1. 找到下一个状态
   。这里proposal distribution设为正态分布。生成
   ,
   .
2. 接受率
   .其中，
   为随机变量
   的概率密度。
3. 接受率的概率实现。如果接受，
   ,否则，
   。

import 

二、Gibbs sampling

1、算法理解

Gibbs sampling适用于高维分布的抽样问题。在M-H抽样算法的基础上，如果我们能够比较容易的得到条件分布，那么就可以通过固定其他维度，一次只对一个维度上的条件分布抽样的方法进行全局抽样。

Gibbs sampling里的接受率恒为1。举例说明，

两个样本点满足马氏链的细致平稳条件。因为

其中，

表示从A点转移到B点的转移概率。所以在二维的分布中，可以得到从任意一个点转移到另外一个点都是平稳的，限制是每次变换只能转移一个维度。二维转移图可如下所示。

2、示例

一只鸡每天会下N个蛋,N服从参数为

的泊松分布，每个鸡蛋成功孵出小鸡的比例为p。p未知，其先验分布服从beta分布。

.参数

已知。我们的观测数据只有每天孵出的小鸡个数

,

属于隐变量，观测不到。如何通过Gibbs Sampling 方法找到p的后验期望呢？

在不引入随机变量N的时候，后验分布比较麻烦。引入N后，可得，

通过迭代，即可得到p,N的抽样值。

x, lambda1, a, b = 7, 10, 1, 1
niter = 10000
p = [0 for i in range(niter)]
N = [0 for i in range(niter)]#初始值
p[0] = 0.5 
N[0] = 2*x
for i in range(1,niter):p[i] = random.betavariate(x+a, N[i-1]-x+b)N[i] = x + np.random.poisson(lambda1*(1-p[i-1]))plt.hist(x = p, bins = 100,normed=True)
plt.hist(x = N, bins = 100,normed=True)  
plt.show()

^[1][2][3]

参考

1. ^https://blog.csdn.net/lin360580306/article/details/51240398
2. ^https://www.jianshu.com/p/dc7624cf6adb
3. ^Introduction to Probability--Joseph K.Blitzstein