MCMC之马尔可夫链
1.马尔可夫链简介
因为某一时刻状态转移只依赖于它的前一个状态，那么我们只要能求出系统中任意两个状态之间的转移概率，进而得到状态转移概率矩阵，那么马尔科夫链的模型便定了。以下图股市模型为例，共有三个状态，分别为牛市(Bull market)、熊市(Bear market)、横盘(Stagnant market)。每一个状态都能够以一定概率转移到下一状态，比如牛市以0.075的概率转移到横盘的概率，这些状态转移概率图可以转换为矩阵的形式进行表示。
如果我们定义矩阵$P$某一位置P(i,j)的值为P(j|i)，即从状态i转移到状态j的概率，并定义牛市的状态为0、熊市状态为1、横盘状态为2，这样便得到马尔可夫链模型的状态转移矩阵。
那么马尔科夫链模型的状态转移矩阵和蒙特卡罗方法所需要的概率分布样本集有什么关系呢？
2.马尔可夫链状态转移矩阵性质
得到马尔可夫链状态转移矩阵，我们看看马尔可夫链模型状态转移矩阵的性质。仍以上面的状态转移矩阵为例，假设当前股市的概率分布为[0.3, 0.4, 0.3]，即30%概率的牛市、40%概率的熊市、30%概率的横盘。如果以这个状态作为序列概率分布的初始状态t0，与状态转移矩阵计算得到t1,t2,t3,…时刻的概率，相应代码和结果如下。
import numpy as np
def markov_chain():
matrix = np.matrix([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05],
[0.25, 0.25, 0.5]], dtype=float)
current = np.matrix([[0.3, 0.4, 0.3]], dtype=float)
for i in range(100):
current = current * matrix
print "Current round:", i + 1
print current
if __name__ == '__main__':
markov_chain()Current round: 1
[[0.405 0.4175 0.1775]]
Current round: 2
[[0.4715 0.40875 0.11975]]
Current round: 3
[[0.5156 0.3923 0.0921]]
Current round: 4
[[0.54591 0.375535 0.078555]]
...
...
Current round: 58
[[0.62499999 0.31250001 0.0625 ]]
Current round: 59
[[0.62499999 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]]
...
...
Current round: 98
[[0.625 0.3125 0.0625]]
Current round: 99
[[0.625 0.3125 0.0625]]
Current round: 100
[[0.625 0.3125 0.0625]]可以发现从第60轮开始，状态转移矩阵的概率分布就不变了，一直保持在[0.625, 0.3125, 0.0625]，即62.5%的概率的牛市、31.25%概率的熊市、6.25%概率的横盘，那么这个是巧合吗？
答案并不是巧合，如果我们换一个初始概率分布t0，比如以[0.7, 0.1, 0.2]作为初始概率分布，然后带入状态转移矩阵得到t1,t2,t3,…时刻的概率，会发现得到同样的结果。即最终的状态概率分布会趋于同一个稳定概率分布[0.625, 0.3125, 0.0625]，也就是说，马尔可夫链的状态转移矩阵收敛到稳定概率分布与初始状态概率分布无关。
上述结果是一个非常好的形式，比如我们得到了稳定概率分布所对应的马尔可夫链模型的状态转移矩阵，那么可以用任意的概率分布样本开始，带入马尔可夫链状态转移矩阵，然后就可以得到符合对应稳定概率分布的样本。这个性质不光对于上面的状态转移矩阵有效，对于绝大多数的其他马尔可夫链模型的状态转移矩阵也有效。同时不光是离散状态，连续状态情况下也成立。
import numpy as np
def markov_chain():
matrix = np.matrix([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05], [0.25, 0.25, 0.5]], dtype=float)
for i in range(10):
matrix = matrix * matrix
print "Current round:", i + 1
print matrix
if __name__ == '__main__':
markov_chain()Current round: 1
[[0.8275 0.13375 0.03875]
[0.2675 0.66375 0.06875]
[0.3875 0.34375 0.26875]]
Current round: 2
[[0.73555 0.212775 0.051675]
[0.42555 0.499975 0.074475]
[0.51675 0.372375 0.110875]]
...
Current round: 5
[[0.62502532 0.31247685 0.06249783]
[0.6249537 0.31254233 0.06250397]
[0.62497828 0.31251986 0.06250186]]
Current round: 6
[[0.625 0.3125 0.0625]
[0.625 0.3125 0.0625]
[0.625 0.3125 0.0625]]
...
Current round: 10
[[0.625 0.3125 0.0625]
[0.625 0.3125 0.0625]
[0.625 0.3125 0.0625]]
上面性质中，有几点需要特意说明
3.基于马尔可夫链采样
4.马尔可夫链总结
如果假定我们可以得到所需要采样样本的平稳分布所对应的马尔可夫链状态转移矩阵，那么我们就可以用马尔可夫链采样得到我们需要的样本集，进而进行蒙特卡罗模拟。但是现在还有个很重要的问题，随意给定一个平稳分布π ，如何得到它所对应的马尔可夫链状态转移矩阵P呢？下篇文章，我们将重点介绍MCMC采用通过与会的方式解决上述问题，以及改进版的M-H采样和Gibbs采样。
你看到的这篇文章来自于公众号「谓之小一」，欢迎关注我阅读更多文章。