敞开生长之旅！这是我参加「日新计划 2 月更文应战」的第 32 天，点击检查活动概况

关于一般的散布的采样，在许多的编程语言中都有完成，如最基本的满意均匀散布的随机数，但是关于复杂的散布，要想对其采样，却没有完成好的函数，在这里，能够使用马尔可夫链蒙特卡罗(Markov Chain Monte Carlo, MCMC)办法，其间Metropolis-Hastings采样和Gibbs采样是MCMC中使用较为广泛的两种形式。

MCMC的基础理论为马尔可夫进程，在MCMC算法中，为了在一个指定的散布上采样，依据马尔可夫进程，首要从任一状况动身，模仿马尔可夫进程，不断进行状况搬运，最终收敛到平稳散布。

1. 马尔可夫链

1.1. 马尔可夫链

设 $X_t$ 表明随机变量 $X$ 在离散时刻 $t$ 时刻的取值。若该变量随时刻变化的搬运概率只是依赖于它的当时取值，即

,Xt=si)=P(Xt+1=sj∣Xt=si)P\left ( X_{t+1}=s_j\mid X_0=s_{0},X_1=s_{1}, \cdots ,X_t=s_i \right )=P\left ( X_{t+1}=s_j\mid X_t=s_i \right )

也就是说状况搬运的概率只依赖于前一个状况。称这个变量为马尔可夫变量，其间， $,si,sj∈s_0,s_1,\cdots ,s_i,s_j\in \Omega$ 为随机变量 $X$ 可能的状况。这个性质称为马尔可夫性质，具有马尔可夫性质的随机进程称为马尔可夫进程。

马尔可夫链指的是在一段时刻内随机变量 $X$ 的取值序列 $,Xm)\left ( X_0,X_1,\cdots ,X_m \right )$ ，它们满意如上的马尔可夫性质。

1.2. 搬运概率

马尔可夫链是经过对应的搬运概率界说的，搬运概率指的是随机变量从一个时刻到下一个时刻，从状况 $s_i$ 搬运到另一个状况 $s_j$ 的概率，即：

P(i→j):=Pi,j=P(Xt+1=sj∣Xt=si)P\left ( i\rightarrow j \right ):=P_{i,j}=P\left ( X_{t+1}=s_j\mid X_t=s_i \right )

记 $k(t)\pi _k^{\left ( t \right )}$ 表明随机变量 $X$ 在时刻 $t$ 的取值为 $s_k$ 的概率，则随机变量 $X$ 在时刻 $t + 1$ 的取值为 $s_i$ 的概率为：

i(t+1)=P(Xt+1=si)=∑kP(Xt+1=si∣Xt=sk)⋅P(Xt=sk)=∑kPk,i⋅k(t)\begin{aligned} \pi _i^{\left ( t+1 \right )} &=P\left ( X_{t+1}=s_i \right ) \\ &= \sum_{k}P\left ( X_{t+1}=s_i\mid X_{t}=s_k \right )\cdot P\left ( X_{t}=s_k \right )\\ &= \sum_{k}P_{k,i}\cdot \pi _k^{\left ( t \right )} \end{aligned}

假设状况的数目为 $n$ ，则有：

,n(t))[P1,1P1,2⋯P1,nP2,1P2,2⋯P2,n⋮⋮⋮Pn,1Pn,2⋯Pn,n]\left ( \pi _1^{\left ( t+1 \right )},\cdots ,\pi _n^{\left ( t+1 \right )} \right )=\left ( \pi _1^{\left ( t \right )},\cdots ,\pi _n^{\left ( t \right )} \right )\begin{bmatrix} P_{1,1} & P_{1,2} & \cdots & P_{1,n}\\ P_{2,1} & P_{2,2} & \cdots & P_{2,n}\\ \vdots & \vdots & & \vdots \\ P_{n,1} & P_{n,2} & \cdots & P_{n,n} \end{bmatrix}

1.3. 马尔可夫链的平稳散布

关于马尔可夫链，需要留意以下的两点：

1、周期性：即经过有限次的状况搬运，又回到了自身；
2、不可约：即两个状况之间相互搬运；

假如一个马尔可夫进程既没有周期性，又不可约，则称为各态遍历的。

关于一个各态遍历的马尔可夫进程，不管初始值 $(0)\pi ^{\left ( 0 \right )}$ 取何值，跟着搬运次数的增多，随机变量的取值散布最终都会收敛到仅有的平稳散布 $∗\pi ^{\ast }$ ，即：

limt→∞(0)Pt=∗\underset{t\rightarrow \infty }{lim}\pi ^{\left ( 0 \right )}\mathbf{P}^t=\pi ^{\ast }

且这个平稳散布 $∗\pi ^{\ast }$ 满意：

∗P=∗\pi ^{\ast }\mathbf{P}=\pi ^{\ast }

其间， $P=(pi,j)nn\mathbf{P}=\left ( p_{i,j} \right )_{n\times n}$ 为搬运概率矩阵。

2. 马尔可夫链蒙特卡罗办法

2.1. 基本思想

关于一个给定的概率散布 $P(X)P\left (X \right )$ ，若是要得到其样本，经过上述的马尔可夫链的概念，我们能够结构一个搬运矩阵为 $P\mathbf{P}$ 的马尔可夫链，使得该马尔可夫链的平稳散布为 $P(X)P\left (X \right )$ ，这样，不管其初始状况为何值，假设记为 $x_0$ ，那么跟着马尔科夫进程的搬运，得到了一系列的状况值，如： $,x_0,x_1,x_2,\cdots ,x_n,x_{n+1},\cdots ,$ ，假如这个马尔可夫进程在第 $n$ 步时已经收敛，那么散布 $P(X)P\left (X \right )$ 的样本即为 $xn,xn+1,⋯x_n,x_{n+1},\cdots$ 。

2.2. 详尽平稳条件

关于一个各态遍历的马尔可夫进程，若其搬运矩阵为 $P\mathbf{P}$ ，散布为 $(x)\pi \left ( x \right )$ ，若满意：

(i)Pi,j=(j)Pj,i\pi \left ( i \right )P_{i,j}=\pi \left ( j \right )P_{j,i}

则 $(x)\pi \left ( x \right )$ 是马尔可夫链的平稳散布，上式称为详尽平稳条件。

2.3. Metropolis采样算法

Metropolis采样算法是最基本的依据MCMC的采样算法。

2.3.1. Metropolis采样算法的基本原理

假设需要从方针概率密度函数 $p()p\left ( \theta \right )$ 中进行采样，一起， $\theta$ 满意 $−∞<<∞-\infty <\theta <\infty$ 。Metropolis采样算法依据马尔可夫链去生成一个序列：

(1)→(2)→⋯(t)→\theta ^{\left ( 1 \right )}\rightarrow \theta ^{\left ( 2 \right )}\rightarrow \cdots \theta ^{\left (t \right )}\rightarrow

其间， $\theta ^{\left (t \right )}$ 表明的是马尔可夫链在第 $t$ 代时的状况。

在Metropolis采样算法的进程中，首要初始化状况值 $(1)\theta ^{\left (1 \right )}$ ，然后利用一个已知的散布 $q(∣(t−1))q\left ( \theta \mid \theta ^{\left ( t-1 \right )} \right )$ 生成一个新的候选状况 $(∗)\theta ^{\left (\ast \right )}$ ，随后依据一定的概率挑选承受这个新值，或者回绝这个新值，在Metropolis采样算法中，概率为：

p((∗))p((t−1)))\alpha =min\: \left ( 1,\; \frac{p\left ( \theta ^{\left ( \ast \right )} \right )}{p\left ( \theta ^{\left ( t-1 \right )} \right )} \right )

这样的进程一向继续到采样进程的收敛，当收敛以后，样本 $(t)\theta ^{\left (t \right )}$ 即为方针散布 $p()p\left ( \theta \right )$ 中的样本。

2.3.2. Metropolis采样算法的流程

依据以上的分析，能够总结出如下的Metropolis采样算法的流程：

初始化时刻 $t = 1$
设置 $u$ 的值，并初始化初始状况 $(t)=u\theta ^{\left (t \right )}=u$
重复一下的进程：
- 令 $t = t + 1$
- 从已知散布 $q(∣(t−1))q\left ( \theta \mid \theta ^{\left ( t-1 \right )} \right )$ 中生成一个候选状况 $(∗)\theta ^{\left (\ast \right )}$
- 计算承受的概率： $p((∗))p((t−1)))\alpha =min\: \left ( 1,\; \frac{p\left ( \theta ^{\left ( \ast \right )} \right )}{p\left ( \theta ^{\left ( t-1 \right )} \right )} \right )$
- 从均匀散布 $Uniform(0,1)Uniform\left ( 0, 1 \right )$ 生成一个随机值 $a$
- 假如 $a⩽a\leqslant \alpha$ ，承受新生成的值： $(t)=(∗)\theta ^{\left (t \right )}=\theta ^{\left (\ast \right )}$ ；否则： $(t)=(t−1)\theta ^{\left (t \right )}=\theta ^{\left (t-1 \right )}$
直到 $t = T$

2.3.3. Metropolis算法的解说

要证明Metropolis采样算法的正确性，最重要的是要证明结构的马尔可夫进程满意如上的详尽平稳条件，即：

(i)Pi,j=(j)Pj,i\pi \left ( i \right )P_{i,j}=\pi \left ( j \right )P_{j,i}

关于上面所述的进程，散布为 $p()p\left ( \theta \right )$ ，从状况 $i$ 搬运到状况 $j$ 的搬运概率为：

Pi,j=i,j⋅Qi,jP_{i,j} =\alpha _{i,j}\cdot Q_{i,j}

其间， $Q_{i,j}$ 为上述已知的散布。

关于挑选该已知的散布，在Metropolis采样算法中，要求该已知的散布有必要是对称的，即 $Q_{i,j}=Q_{j,i}$ ，即 $q(=(t)∣(t−1))=q(=(t−1)∣(t))q\left ( \theta =\theta ^{\left ( t \right )}\mid \theta ^{\left ( t-1 \right )} \right )=q\left ( \theta =\theta ^{\left ( t-1 \right )}\mid \theta ^{\left ( t \right )} \right )$ 常用的符合对称的散布主要有：正态散布，柯西散布以及均匀散布等。

接下来，需要证明在Metropolis采样算法中结构的马尔可夫链满意详尽平稳条件。

{p((i))p((j)),1}⋅Qj,i=p((j))⋅j,i⋅Qj,i=p((j))Pj,i\begin{aligned} p\left ( \theta ^{\left ( i \right )} \right )P_{i,j} &=p\left ( \theta ^{\left ( i \right )} \right )\cdot \alpha _{i,j}\cdot Q_{i,j} \\ &= p\left ( \theta ^{\left ( i \right )} \right )\cdot min\; \left \{ 1,\frac{p\left ( \theta ^{\left ( j \right )} \right )}{p\left ( \theta ^{\left ( i \right )} \right )} \right \}\cdot Q_{i,j}\\ &=min\; \left \{ p\left ( \theta ^{\left ( i \right )} \right )Q_{i,j},p\left ( \theta ^{\left ( j \right )} \right )Q_{i,j} \right \}\\ &=p\left ( \theta ^{\left ( j \right )} \right )\cdot min\; \left \{ \frac{p\left ( \theta ^{\left ( i \right )} \right )}{p\left ( \theta ^{\left ( j \right )} \right )}, 1 \right \}\cdot Q_{j,i}\\ &=p\left ( \theta ^{\left ( j \right )} \right )\cdot \alpha _{j,i}\cdot Q_{j,i}\\ &=p\left ( \theta ^{\left ( j \right )} \right )P_{j,i} \end{aligned}

因而，经过以上的办法结构出来的马尔可夫链是满意详尽平稳条件的。

2.3.4. 试验

假设需要从柯西散布中采样数据，我们利用Metropolis采样算法来生成样本，其间，柯西散布的概率密度函数为：

f()=1(1+2)f\left ( \theta \right )=\frac{1}{\pi \left ( 1+\theta ^2 \right )}

那么，依据上述的Metropolis采样算法的流程，承受概率 $\alpha$ 的值为：

(1,1+[(t)]21+[(∗)]2)\alpha =min\; \left ( 1,\frac{1+\left [ \theta ^{\left ( t \right )} \right ]^2}{1+\left [ \theta ^{\left ( \ast \right )} \right ]^2} \right )

代码如下：

'''
Date:20160629
@author: zhaozhiyong
'''
import random
from scipy.stats import norm
import matplotlib.pyplot as plt
def cauchy(theta):
    y = 1.0 / (1.0 + theta ** 2)
    return y
T = 5000
sigma = 1
thetamin = -30
thetamax = 30
theta = [0.0] * (T+1)
theta[0] = random.uniform(thetamin, thetamax)
t = 0
while t < T:
    t = t + 1
    theta_star = norm.rvs(loc=theta[t - 1], scale=sigma, size=1, random_state=None)
    #print theta_star
    alpha = min(1, (cauchy(theta_star[0]) / cauchy(theta[t - 1])))
    u = random.uniform(0, 1)
    if u <= alpha:
        theta[t] = theta_star[0]
    else:
        theta[t] = theta[t - 1]
ax1 = plt.subplot(211)
ax2 = plt.subplot(212) 
plt.sca(ax1)
plt.ylim(thetamin, thetamax)
plt.plot(range(T+1), theta, 'g-')
plt.sca(ax2)
num_bins = 50
plt.hist(theta, num_bins, normed=1, facecolor='red', alpha=0.5)
plt.show()

试验的成果：

关于Metropolis采样算法，其要求选定的散布有必要是对称的，为了补偿这样的一个缺点，鄙人一篇中，介绍一下Metropolis-Hastings采样算法，其是Metropolis采样算法的推广形式。

参考文献

[1] 马尔可夫链蒙特卡罗算法

[2] 受限玻尔兹曼机（RBM）学习笔记（一）准备常识

[3] LDA数学八卦

敞开生长之旅！这是我参加「日新计划 2 月更文应战」的第 32 天，点击检查活动概况

马尔可夫链蒙特卡罗方法MCMC

1. 马尔可夫链

1.1. 马尔可夫链

1.2. 搬运概率

1.3. 马尔可夫链的平稳散布

2. 马尔可夫链蒙特卡罗办法

2.1. 基本思想

2.2. 详尽平稳条件

2.3. Metropolis采样算法

2.3.1. Metropolis采样算法的基本原理

2.3.2. Metropolis采样算法的流程

2.3.3. Metropolis算法的解说

2.3.4. 试验

参考文献

作者信息

推广

马尔可夫链蒙特卡罗方法MCMC

1. 马尔可夫链

1.1. 马尔可夫链

1.2. 搬运概率

1.3. 马尔可夫链的平稳散布

2. 马尔可夫链蒙特卡罗办法

2.1. 基本思想

2.2. 详尽平稳条件

2.3. Metropolis采样算法

2.3.1. Metropolis采样算法的基本原理

2.3.2. Metropolis采样算法的流程

2.3.3. Metropolis算法的解说

2.3.4. 试验

参考文献

相关文章

异步超时中断，知其然，也要知其所以然～

图片处理看这篇就完了「GitHub 热点速览 v.21.48」

到底应不应该使用 lombok

羊羊刷题笔记Day13/60 | 第六章 二叉树P1 | 二叉树理论基础、递归遍历、迭代遍历与总结

作者信息

推广

羊羊刷题笔记Day13/60 | 第六章二叉树P1 | 二叉树理论基础、递归遍历、迭代遍历与总结