马尔科夫链及其采样方法

前言

了解了什么是蒙特卡罗方法之后，自然引出了马尔可夫链这个概念，其在 WIKI 的解释如下：

马尔可夫链（英语：Markov chain），又称离散时间马尔可夫链（discrete-time Markov chain，缩写为DTMC），因俄国数学家安德烈·马尔可夫得名，为状态空间中经过从一个状态到另一个状态的转换的随机过程。该过程要求具备“无记忆”的性质：下一状态的概率分布只能由当前状态决定，在时间序列中它前面的事件均与之无关。这种特定类型的“无记忆性”称作马尔可夫性质。 在马尔可夫链的每一步，系统根据概率分布，可以从一个状态变到另一个状态，也可以保持当前状态。状态的改变叫做转移，与不同的状态改变相关的概率叫做转移概率。

定义及性质

依据上述 wiki 的定义可以知道，马尔可夫链能在很多时间序列模型中得到广泛的应用，比如循环神经网络RNN，隐式马尔科夫模型HMM等，当然MCMC也需要它。

如果用精确的数学定义来描述，则假设我们的序列状态是$...X_{t-2}, X_{t-1}, X_{t}, X_{t+1},...$那么我们的在时刻 $X_{t+1}$ 的状态的条件概率仅仅依赖于时刻 $X_t$，即：

$$P(X_{t+1} |...X_{t-2}, X_{t-1}, X_{t} ) = P(X_{t+1} | X_{t})$$

既然某一时刻状态转移的概率只依赖于它的前一个状态，那么我们只要能求出系统中任意两个状态之间的转换概率，这个马尔科夫链的模型就定了。我们来看看下图这个马尔科夫链模型的具体的例子(来源于维基百科)。

这个马尔科夫链是表示股市模型的，共有三种状态：牛市（Bull market）, 熊市（Bear market）和横盘（Stagnant market）。每一个状态都以一定的概率转化到下一个状态。比如，牛市以0.025的概率转化到横盘的状态。这个状态概率转化图可以以矩阵的形式表示。如果我们定义矩阵P某一位置$P(i,j)$的值为$P(j|i)$,即从状态i转化到状态j的概率，并定义牛市为状态0， 熊市为状态1, 横盘为状态2. 这样我们得到了马尔科夫链模型的状态转移矩阵为：

$$P=\left( \begin{array}{ccc} 0.9&0.075&0.025 \\ 0.15&0.80& 0.05 \\ 0.25&0.25&0.5 \end{array} \right)$$

假设已知当前处在熊市、牛市、横盘的人的比例是概率分布向量$π_0=[π_0(1),π_0(2),π_0(3)]$,那么第一天的分布比例将是$π_1=π_0P$, 第二天的分布比例将是$π_2=π_1P=π_0P^2$ ，第n天分布比例将是$π_n=π_{n−1}P=π_0P^n$。

假设我们当前已知股市的概率分布为：[0.8,0.1,0.1],即80%概率的牛市，10%概率的熊盘与10%的横盘。然后这个状态作为序列概率分布的初始状态$t_0$，将其带入这个状态转移矩阵计算$t_1,t_2,t_3...$的状态。代码如下：

matrix = np.matrix([[0.9,0.075,0.025],[0.15,0.8,0.05],[0.25,0.25,0.5]], dtype=float)
vector1 = np.matrix([[0.8,0.1,0.1]], dtype=float)
for i in range(100):
    vector1 = vector1*matrix
    print "Current round:" , i+1
    print vector1

通过观察输出结果,可以看出在迭代一定次数之后,数据结果出现了如下的稳定情况

Current round: 48
[[0.62500011 0.3124999  0.06249999]]
Current round: 49
[[0.62500008 0.31249992 0.06249999]]
Current round: 50
[[0.62500006 0.31249994 0.06249999]]
.............
Current round: 59
[[0.625  0.3125 0.0625]]
Current round: 60
[[0.625  0.3125 0.0625]]
Current round: 61
[[0.625  0.3125 0.0625]]
Current round: 62
[[0.625  0.3125 0.0625]]
........

可以发现，从第60轮开始，我们的状态概率分布就不变了，一直保持在[0.625 0.3125 0.0625]，即62.5%的牛市，31.25%的熊市与6.25%的横盘。

我们现在换一个初始概率分布试一试，现在我们用[0.7,0.1,0.2]作为初始概率分布, 利用上述的代码计算,得出的结果一样出现的最终状态的概率分布趋于同一个稳定的概率分布[0.625 0.3125 0.0625], 也就是说我们的马尔科夫链模型的状态转移矩阵收敛到的稳定概率分布与我们的初始状态概率分布无关。这是一个非常好的性质, 也就是说，如果我们得到了这个稳定概率分布对应的马尔科夫链模型的状态转移矩阵，则我们可以用任意的概率分布样本开始，代入马尔科夫链模型的状态转移矩阵，这样经过一些序列的转换，最终就可以得到符合对应稳定概率分布的样本。

这个性质不光对我们上面的状态转移矩阵有效，对于绝大多数的其他的马尔科夫链模型的状态转移矩阵也有效。同时不光是离散状态，连续状态时也成立。

同时这个性质也意味着，对于一个确定的状态转移矩阵 P，它的 n 次幂$P^n$在当 n 大于一定的值的时候也可以发现是确定的.

性质的数学语言描述

如果一个非周期的马尔科夫链有状态转移矩阵P, 并且它的任何两个状态是连通的，那么$\lim_{n \to \infty}P_{ij}^n$与 i 无关，我们有：

$$\lim_{n \to \infty}P_{ij}^n = \pi(j) \tag{1}$$ $$\lim_{n \to \infty}P^n = \left( \begin{array}{ccc} \pi(1)&\pi(2)&\ldots&\pi(j)&\ldots \\ \pi(1)&\pi(2)&\ldots&\pi(j)&\ldots \\ \ldots&\ldots&\ldots&\ldots&\ldots \\ \pi(1)&\pi(2)&\ldots&\pi(j)&\ldots \\ \ldots&\ldots&\ldots&\ldots&\ldots \end{array} \right) \tag{2}$$ $$\pi(j) = \sum\limits_{i=0}^{\infty}\pi(i)P_{ij} \tag{3}$$

π 是方程$πP=π$的唯一非负解，其中

$$\pi = [\pi(1),\pi(2),...,\pi(j),...]\;\; , \sum\limits_{i=0}^{\infty}\pi(i) = 1 \tag{4}$$

上面的性质中需要解释的有：

1. 非周期的马尔科夫链：定理中的“非周期“这个概念不解释，因为我们遇到的绝大多数马氏链都是非周期的；这个主要是指马尔科夫链的状态转化不是循环的，如果是循环的则永远不会收敛。幸运的是我们遇到的马尔科夫链一般都是非周期性的。用数学方式表述则是：对于任意某一状态i，d为集合$\{n \mid n \geq 1,P_{ii}^n>0 \}$的最大公约数，如果 d=1 ，则该状态为非周期的
2. 任何两个状态是连通的：这个指的是从任意一个状态可以通过有限步到达其他的任意一个状态，不会出现条件概率一直为0导致不可达的情况。
3. 马尔科夫链的状态数可以是有限的，也可以是无限的。因此可以用于连续概率分布和离散概率分布。
4. π 通常称为马尔科夫链的平稳分布。
5. 我们用$X_i$表示在马氏链上跳转第 i 步所处的状态，如果$\lim_{n \rightarrow \infty}P_{ij}^{n} = \pi(j)$存在,很容易证明定理第二个结论

$$P(X_{n+1}=j) = \sum_{i=1}^{\infty}P(X_{n}=i)P(X_{n+1}=j\mid x_{n}=i)=\sum_{i=1}^{\infty}P(X_{n}=i) p_{ij}$$

基于马尔科夫链采样

对于给定的概率分布$p(x)$ ，我们希望能有便捷的方式生成它对应的样本。由于马尔科夫链能收敛到平稳分布，于是一个很的漂亮想法是：如果我们能构造一个转移矩阵为 P 的马尔科夫链，使得该马尔科夫链的平稳分布恰好是$p(x)$，那么我们从任何一个初始状态$x_0$出发沿着马尔科夫链转移，得到一个转移序列$x_0,x_1,x_2,\ldots,x_n,x_{n+1},\ldots$, 如果马尔科夫在第n步已经收敛了，于是我们就得到了$p(x)$的样本$x_n,x_{n+1},\ldots$

总结下基于马尔科夫链的采样过程：

1. 输入马尔科夫链状态转移矩阵P，设定状态转移次数阈值$n_1$，需要的样本个数$n_2$
2. 从任意简单概率分布采样得到初始状态值$x_0$
3. $for t=0 to n_1+n_2−1:$ 从条件概率分布$P(x|x_t)$中采样得到样本$x_t+1$

样本集$(x_{n1},x_{n_1+1},\ldots,...,x_{n1+n2−1})$即为我们需要的平稳分布对应的样本集。

小结

马尔科夫链的收敛性质主要由转移矩阵P决定， 所以基于马尔科夫链做采样的关键问题是如何构造转移矩阵P，使得平稳分布恰好是我们要的分布$p(x)$。

如何能做到这一点呢？

幸运的是，MCMC采样通过迂回的方式解决了上面这个大问题，我们在下一篇来讨论MCMC的采样，以及它的使用改进版采样: M-H(Metropolis-Hastings)采样和Gibbs采样.