蒙特卡罗方法和马尔科夫链简明文档

初次编辑于2020-06-09，加入蒙特卡罗方法以及马尔科夫矩阵的介绍

蒙特卡罗方法

蒙特卡罗方法是一种随机的采样方法，所求解的问题同一定的概率模型相联系，用电子计算机实现统计模拟或以获得问题的近似解。

两种应用
蒙特卡罗方法有两种应用：一种是计算复杂图形的面积；一种是计算复杂函数的积分。但是从原理上来说，其实都是再求复杂图形的面积。

计算圆周率.gif

假设我们有一袋豆子，均匀的撒在如上图的正方形区域内，通过数出落在红色区域内的豆子数和总的豆子数作比，我们可以近似得求得圆周率 $\pi$ ；并且撒出去的豆子越多，求得圆周率越接近真实值。

计算复杂定积分.png

蒙特卡罗方法的另外一个应用就是计算上图所示函数 $f(x)$ 在闭区间 $[a, b]$ 的积分 $\int_a^bf(x)dx$ ，亦即计算图中 $f(x)$ 和坐标轴所围阴影面积。
计算上述定积分一个比较简单粗暴的方法是任取一 $x_0\in[a, b]$ ，定积分的值为：

$(b - a)f(x_0)$

当然，用一个值代表[a,b]区间上所有的𝑓(𝑥)的值，这太过于简单粗暴了。那么我们可以采样 $[a,b]$ 区间的n个值： $𝑥_0,𝑥_1,...𝑥_{𝑛−1}$ ,用它们的均值来代表 $[a,b]$ 区间上所有的 $𝑓(𝑥)$ 的值。这样，上面的定积分的近似求解为:

$\frac{b-a}{n}\sum_{i=0}^{n-1}f(x_i)$

虽然上面的方法可以一定程度上求解出近似的解，但是它隐含了一个假定，即𝑥在[a,b]之间是均匀分布的，而绝大部分情况，𝑥在[a,b]之间不是均匀分布的。如果采用上面的方法，则模拟求出的结果很可能和真实值相差甚远。

为了解决这个问题，可以假定得到x在[a,b]之间的概率分布 $p(x)$ ，这样我们的定积分求和如下所示：

$\theta = \int_a^bf(x)dx = \int_a^b\frac{f(x)}{p(x)}p(x)dx \approx \frac{1}{n}\sum_{i=0}^{n-1}\frac{f(x_i)}{p(x_i)}$

上式的最后一项，相当于离散值求分布 $\frac{f(x)}{p(x)}$ 期望；倒数第二项，相当于求连续值分布 $\frac{f(x)}{p(x)}$ 的期望。上式最右边的这个形式就是蒙特卡罗方法的一般形式。当然这里是连续函数形式的蒙特卡罗方法，但是在离散时一样成立。

至此，我们只需要知道随机密度函数 $p(x)$ 就可以通过蒙特卡洛方法近似计算发杂积分值。而借助马尔科夫链，我们无需知道具体的先验概率分布，只需要知道状态转换的概率矩阵，我们就可以求得后验概率分布，从而创造使用蒙特卡罗方法的条件。

马尔科夫链

定义

马尔科夫链的定义是，假设某一时刻状态转移的概率只依赖于它的前一刻状态。比如假设每天的天气是一个状态，那么今天是不是晴天只依赖于昨天的天气，而与前天甚至更早前的天气没有任何关系。这样假设比较草率，但是这样做可以大大简化模型的复杂度。

如果要用精确的数学语言来描述，那么假设一组序列状态为 $...X_{t-2},X_{t_1}, X_t, X_{t+1}, X_{t+2}...$ 根据马尔科夫链的定义，状态 $X_{t+1}$ 将只依赖于其前一刻的状态 $X_t$ ，即
$P(X_{T+1}|...X_{t-2}, X_{t_1}, X_t) = P(X_{T+1}|X_t)$

由于某一时刻的状态转移概率只依赖于其前一刻的状态，所以只要确定系统中任意两个状态之间的转换概率，这个马尔科夫链的模型也就确定了。下面以一个股市状态例子来进行说明。

1024px-Financial_Markov_process.svg.png

说明：图中的箭头是状态之间的转换；数字表明状态之间转换的概率；Bull market、Bear Market以及Stagnant Market依次表示牛市、熊市和横盘。

每一个状态都以一定的概率转换到下一个状态，比如牛市有2%的概率转换到熊市。由此，这个状态概率转换图可以以矩阵的形式进行标示。假设以 $P$ 来表示这个矩阵，位置 $P(i, j)$ 的值为 $P(i|j)$ 表示从状态j转换到状态i的概率，并且假设牛市（Bull Market）为状态0，熊市（Bear market）为状态1，横盘（Stagnant Market）为状态0，那么对应该图的状态转换矩阵 $P$ 如下所示：
$P = \begin{pmatrix} 0.9 & 0.075 & 0.025 \\\\ 0.15 & 0.8 & 0.05 \\\\ 0.25 & 0.25 & 0.5 \end{pmatrix}$

性质

下面通过几个例子来看一下马尔科夫链的性质。
假设当前股市的概率为：[0.3, 0.4, 0.3]，即当前股市有30%的概率为牛市，40%的概率为熊市，40%的概率横盘。以此状态为 $t_0$ 时刻的初始状态，然后利用状态转换矩阵 $P$ 计算 $t_1,t_2,t_3...$ 时刻的状态。代码如下：

import numpy as np

P =  np.array([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05], [0.25, 0.25, 0.5]], dtype=np.float32)
state = np.array([[0.3, 0.4, 0.3]], dtype=np.float32)

for i in range(60):
  state = np.dot(state, p)
  print("Time : ", i)
  print(state)

部分的输出结果如下：

Time step:  0
     [[0.40500003 0.41750002 0.17750001]]
Time step:  1
     [[0.47150004 0.40875003 0.11975001]]
Time step:  2
     [[0.5156     0.39230004 0.09210001]]
Time step:  3
     [[0.54591    0.37553504 0.078555  ]]
Time step:  4
     [[0.56728804 0.36101004 0.071702  ]]
Time step:  5
     [[0.58263624 0.34928015 0.0680837 ]]
Time step:  6
     [[0.5937855  0.3401428  0.06607176]]
Time step:  7
     [[0.6019463  0.3331661  0.06488766]]
Time step:  8
     [[0.6079485  0.32790074 0.0641508 ]]
Time step:  9
     [[0.6123764  0.32395446 0.06366915]]
Time step:  10
     [[0.6156492  0.3210091  0.06334171]]
...
Time step:  43
     [[0.62499946 0.31250048 0.06250004]]
Time step:  44
     [[0.6249996  0.31250036 0.06250003]]
Time step:  45
     [[0.62499964 0.31250027 0.06250002]]
Time step:  46
     [[0.6249997  0.31250018 0.06250001]]
Time step:  47
     [[0.62499976 0.31250012 0.06250001]]
Time step:  48
     [[0.62499976 0.31250006 0.0625    ]]
Time step:  49
     [[0.62499976 0.31250006 0.0625    ]]
Time step:  50
     [[0.62499976 0.31250006 0.0625    ]]

可以看到，从第48个时间步开始，状态转换概率矩阵就稳定在了 $[0.625, 0.3125, 0.0625]$ ，即转换至下一时刻状态时，有62.5%的概率转换为牛市，有31.25%的概率转换为熊市，有6.25%的机会转换为横盘。

现在我们将初始时刻 $t_0$ 的状态转换概率设置为 $[0.6, 0.1, 0.3]$ ，再次运行前述代码段，得到 $t_1,t_2,t_3...$ 时刻的状态转换概率，部分输出结果如下所示：

Time step:  0
     [[0.63 0.2  0.17]]
Time step:  1
     [[0.6395     0.24975002 0.11075   ]]
Time step:  2
     [[0.64070004 0.27545002 0.08385   ]]
Time step:  3
     [[0.63891    0.28937504 0.071715  ]]
Time step:  4
     [[0.636354   0.29734704 0.06629901]]
Time step:  5
     [[0.6338954  0.30217892 0.06392571]]
Time step:  6
     [[0.6318141  0.30526674 0.06291918]]
Time step:  7
     [[0.6301525  0.30732924 0.06251828]]
Time step:  8
     [[0.6288662  0.30875438 0.06237942]]
Time step:  9
     [[0.62788755 0.3097633  0.06234908]]
Time step:  10
     [[0.62715054 0.31048948 0.0623599 ]]
...
Time step:  40
     [[0.625      0.31249964 0.06249995]]
Time step:  41
     [[0.62499994 0.3124997  0.06249996]]
Time step:  42
     [[0.6249999  0.31249976 0.06249996]]
Time step:  43
     [[0.6249998  0.3124998  0.06249997]]
Time step:  44
     [[0.62499976 0.31249982 0.06249997]]
Time step:  45
     [[0.62499976 0.31249985 0.06249997]]
Time step:  46
     [[0.62499976 0.31249988 0.06249997]]
Time step:  47
     [[0.62499976 0.31249988 0.06249997]]
Time step:  48
     [[0.62499976 0.31249988 0.06249998]]
Time step:  49
     [[0.62499976 0.31249988 0.06249998]]
Time step:  50
     [[0.62499976 0.31249988 0.06249998]]

可以看到，从第44个时间步开始，状态转换概率矩阵就稳定在了 $[0.625， 0.3125， 0.0625]$ 。两次不同的初始 $t_0$ 时刻的状态转换概率矩阵，经过有限次的时间步之后，都会收敛到一个固定的概率转换矩阵上，也就是说马尔科夫链模型的状态转移矩阵收敛到的稳定概率分布与我们的初始状态概率分布无关。这是一个非常好的性质，也就是说，如果得到了这个稳定概率分布对应的马尔科夫链模型的状态转移矩阵，则可以用任意的概率分布样本开始，带入马尔科夫链模型的状态转移矩阵，这样经过一些序列的转换，最终就可以得到符合对应稳定概率分布的样本。

同时，对于一个确定的状态转移矩阵𝑃，它的n次幂 $P^n$ 在当n大于一定的值的时候也可以发现是确定的，以前述矩阵 $P$ 为例，示例计算代码如下所示：

import numpy as np

p = np.array([[0.9, 0.075, 0.025], [0.15, 0.8, 0.05], [0.25, 0.25, 0.5]], dtype=np.float32)
t_p = p 
for i in range(60):
  t_p = np.dot(t_p, p)
  print("Time step: ", i)
  print(t_p)

部分输出如下所示：

Time step:  0
[[0.8275  0.13375 0.03875]
 [0.2675  0.66375 0.06875]
 [0.3875  0.34375 0.26875]]
Time step:  1
[[0.77449995 0.17875001 0.04675   ]
 [0.35750002 0.56825    0.07425   ]
 [0.46749997 0.37125003 0.16125001]]
Time step:  2
[[0.7355499  0.212775   0.051675  ]
 [0.42555    0.499975   0.07447501]
 [0.51675    0.37237504 0.11087501]]
Time step:  3
[[0.70682997 0.23830502 0.054865  ]
 [0.47661003 0.45051503 0.07287501]
 [0.54865    0.36437505 0.08697501]]
Time step:  4
[[0.6856089  0.2573725  0.0570185 ]
 [0.514745   0.41437653 0.07087851]
 [0.57018507 0.35439256 0.07542251]]
Time step:  5
[[0.6699085  0.2715733  0.0585181 ]
 [0.5431466  0.38782674 0.06902671]
 [0.58518106 0.34513357 0.06968551]]
...
Time step:  45
[[0.6250001  0.3124997  0.06249996]
 [0.6249993  0.31250042 0.06250003]
 [0.6249998  0.31250036 0.06250004]]
Time step:  46
[[0.62500006 0.31249976 0.06249997]
 [0.6249994  0.3125003  0.06250001]
 [0.6249999  0.31250027 0.06250004]]
Time step:  47
[[0.625      0.31249982 0.06249998]
 [0.6249995  0.3125002  0.06250001]
 [0.62499994 0.3125002  0.06250003]]
Time step:  48
[[0.625      0.31249985 0.06249998]
 [0.62499964 0.31250015 0.0625    ]
 [0.625      0.31250018 0.06250003]]
Time step:  49
[[0.625      0.31249988 0.06249998]
 [0.62499964 0.3125001  0.0625    ]
 [0.625      0.31250015 0.06250002]]
Time step:  50
[[0.625      0.3124999  0.06249999]
 [0.62499964 0.31250006 0.0625    ]
 [0.625      0.31250012 0.06250001]]
Time step:  51
[[0.625      0.31249994 0.06249999]
 [0.62499964 0.31250003 0.06249999]
 [0.625      0.3125001  0.06250001]]

可以看到，当 $n\ge48$ 的时候， $P^n$ 的值稳定下来，不再变化，并且每一行都是 $[0.625, 0.3125, 0.0625]$ ，这与前述的稳定分布是一致的。据此，可以总结马尔科夫链的收敛性质了。