【翻译活动】面向数据科学的概率论-10.马尔科夫链

原文：prob140/textbook/notebooks/ch_10

译者：喵十八

协议：CC BY-NC-SA 4.0

自豪地采用谷歌翻译

术语说明

条件概率分布（Conditional Probability Distribution，或者条件分布，Conditional Distribution）是现代概率论中的概念：已知两个相关的随机变量X和Y，随机变量Y在条件{X=x}下的条件概率分布是指当已知X的取值为某个特定值x之时，Y的概率分布。 如果Y在条件{X=x}下的条件概率分布是连续分布，那么其密度函数称作Y在条件{X=x}下的条件概率密度函数（条件分布密度、条件密度函数）。与条件分布有关的概念，常常以“条件”作为前缀，如条件期望、条件方差等等。

转移 与转移概率：从状态1变为状态2，称之为状态转移，其对应的概率称之为转移概率。

可数无穷：是指集合中的元素可以与自然数一一对应,也就是说可以用自然数来”数”它的数量,从而其数量为可数无穷.

本章所需python包

# HIDDEN
from datascience import *
from prob140 import *
import numpy as np
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
%matplotlib inline
import math
from scipy import stats
from scipy import misc

# HIDDEN
s = np.arange(1, 6)

def refl_walk_probs(i, j):
    # staying in the same state
    if i-j == 0:
        return 0.5
    
    # moving left or right
    elif 2 <= i <= 4:
        if abs(i-j) == 1:
            return 0.25
        else:
            return 0
        
    # moving right from 1
    elif i == 1:
        if j == 2:
            return 0.5
        else:
            return 0
    
    # moving left from 5
    elif i == 5:
        if j == 4:
            return 0.5
        else:
            return 0

reflecting_walk = MarkovChain.from_transition_function(s, refl_walk_probs)

马尔科夫链

一个 随机过程 是某一概率空间（这里的空间，理解为域更合适，是指一系列随机状态的合集。不是与时间相对的空间，后面的时间也类似）中一系列随机状态的集合。我们将研究一种在离散时间域内演化过程，即有如下随机状态 $X_0, X_1, X_2, \ldots $。想象一下，从时刻0的状态 $X_0$ 开始,不断执行对应时间的操作，状态随之转移。以此类推，在时刻$n$ 时，转移为状态 $n$。

我们已经见过这种过程的例子。例如伯努利实验中，一系列的抛硬币事件构成的独立同分布的序列，就形成这样的过程。每次事件在0和1两个值之间来回传递，每个值都独立于所有其他值。但在许多有趣的过程中，未来的事件的值依赖于当前事件的值，以及过去事件的值。我们可以使用过去和现在来预测未来。

马尔科夫链，是以安德烈·马尔科夫的名字命名的一类随机过程。一个非正式的描述如下，对于马尔科夫链而言，将来的状态的值只取决于当前状态的值，而与如何达到当前状态的值无关。这被称为马尔科夫性质。从形式上看

• 对于任一$n \ge 1$， $X_{n+1}$ 的条件分布，只取决于 $X_0, X_1, \ldots , X_n$中的$X_n$
• 也就是说，对于每一个可能值的序列$i0, i_1, \ldots, i_n, i{n+1}$,

$$P(X_{n+1} = i_{n+1} \mid X_0 = i_0, X_1 = i_1 , \ldots, X_{n-1} = i_{n-1}, X_n = i_n) = P(X_{n+1} = i_{n+1} \mid X_n = i_n)$$

例如在一个随机漫步试验中，赌徒以a美元的财富开始，然后连续投掷一枚公平的硬币（正反面概率都为50%）。如果硬币为正面，他就获得1美元，如果是反面，他就输掉1美元。
设$X{0} = a$，则$n > 0$令$X{n+1} = X_n + I_n$，其中$I_1, I_2, \ldots $是伯努利实验中的独立同分布序列。马尔科夫性质适用于整个过程：给出赌徒在时刻$n$的财富数，那他在时刻$n+1$时的财富数的取值只与其在时刻$n$时的财富数有关，与$n$之前无关。所以该过程$X_0, X_1, X_2, \ldots $是一个马尔科夫链，代表着赌徒的财富随时间的演变。

马尔可夫链的状态空间是链中状态可能值的集合。上述随机漫步的状态空间是所有整数的集合。在本课程中，我们将状态空间限制为离散且有限的。

条件独立

两个随机变量$X$和$Y$是相互独立是指，$X$处于条件$Y$的情况的条件分布和不处于$Y$的情况下的条件分布是一致的。
随机变量$X$和$Y$相对于$Z$条件独立 是指，$X$在条件$Y$和$Z$的情况下的条件分布和在条件$Z$的情况下的条件分布式一致。也就是说，$Z$这一关于$Y$的额外条件，不会影响$X$的值。

在马尔可夫链中，如果定义时刻$n$为现在，定义时刻$n+1$为未来，时刻序列$0$到$n-1$作为过去，马尔科夫性质意味着过去和未来是条件独立的。

转移

设$X_0, X_1, X_2, \ldots $为状态空间$S$中马尔科夫链。根据马尔科夫性质，有限长度的路径或轨迹的概率如下：

$$\begin{align*} & P(X_0 = i_0, X_1 = i_1, X_2 = i_2, \ldots, X_n = i_n) \\ & = ~ P(X_0 = i_0)P(X_1 = i_1 \mid X_0 = i_0)P(X_2 = i_2 \mid X_1 = i_1) \cdots P(X_n = i_n \mid X_{n-1} = i_{n-1}) \end{align*}$$

上式中的条件概率，也被称为转移概率。对于状态$i$和$j$，条件概率$P(X_{n+1} = j \mid X_n = i)$被称为时刻$n$时的一阶转移概率。

对于许多链，例如随机漫步，这些一阶转移概率仅仅由状态$i$和$j$决定，而与时刻$n$无关。
示例：

$$\begin{equation} P(X_{n+1} = j \mid X_n = i) = \begin{cases} \frac{1}{2} & \text{if } j = i-1 \text{ or } j = i+1 \\ 0 & \text{ otherwise} \end{cases} \end{equation}$$

对所有的$n$都与时刻无关。

固定转移概率

当一阶转移概率与时刻$n$无关时，称之为固定或者时间同质的。我们将在本课程中学习的所有马尔可夫链都具有时间同质的转移概率。
对于这样的链，定义一阶转移概率如下：

$$P(i, j) ~ = ~ P(X_{n+1} = j \mid X_n = i) ~ = ~ P(X_1 = j \mid X_0 = i)$$

Then

$$P(X_0 = i_0, X_1 = i_1, X_2 = i_2, \ldots, X_n = i_n) ~ = ~ P(X_0 = i_0)P(i_0, i_1)P(i_1, i_2) \cdots P(i_{n-1}, i_n)$$

一阶转移概率可以表示为矩阵的元素。这不仅仅是为了符号的紧凑-它导致了一个强大的理论。

一阶转移矩阵

链的一阶转移矩阵描述如下，矩阵$\mathbb{P}$，其中$(i, j)$位置处的元素是$P(i, j) = P(X_1 = j \mid X_0 = i)$。

通常，$\mathbb{P}$简称为转移矩阵。注意两个重要属性：

• $\mathbb{P}$是一个正方形矩阵: 它的行和列都由状态空间索引构成。
• $\mathbb{P}$的每一行: 对任一状态$i$, 和时刻$n$, 行$i$ 包含了在$Xn = i$ 情况下，$X{n+1}$的条件分布。 因为它的每一行的和都为1， $\mathbb{P}$ 也被称为 随机矩阵.

让我们看一下示例中转移矩阵的样子。

粘性反转随机漫步

通常，马尔可夫链的转移行为更容易在转移图而不是矩阵中描述。下面是状态1,2,3,4和5上的链的转移图。该图显示了一阶转移概率。

• 如果链条处于任何状态，它移动到原有状态的概率为0.5。
• 如果链处于状态2到4，则它移动到其两个相邻状态中的一个的概率为0.25。
• 如果链处于状态1或5，则它移动到其相邻状态的概率为0.5。

我们称其为反转是在状态1和5可以反转掉头进行转移。整个漫步过程有粘性是指其可能移动到原有状态。

转移图非常适合理解链移动的规则。但是，对于计算，转移矩阵更有帮助。

要开始构造矩阵，我们将数组s设置为状态集，并为转移函数refl_walk_probs 设置成入参为$i$和$j$，返回值为$P(i, j)$的形式。

s = np.arange(1, 6)

def refl_walk_probs(i, j):
    # staying in the same state
    if i-j == 0:
        return 0.5
    
    # moving left or right
    elif 2 <= i <= 4:
        if abs(i-j) == 1:
            return 0.25
        else:
            return 0
        
    # moving right from 1
    elif i == 1:
        if j == 2:
            return 0.5
        else:
            return 0
    
    # moving left from 5
    elif i == 5:
        if j == 4:
            return 0.5
        else:
            return 0

您可以使用prob140库来构造MarkovChain对象。from_transition_function方法有两个参数：

• 状态构成的数组
• 转移函数

并显示MarkovChain对象的一阶转移矩阵。

reflecting_walk = MarkovChain.from_transition_function(s, refl_walk_probs)
reflecting_walk

	1	2	3	4	5
1	0.50	0.50	0.00	0.00	0.00
2	0.25	0.50	0.25	0.00	0.00
3	0.00	0.25	0.50	0.25	0.00
4	0.00	0.00	0.25	0.50	0.25
5	0.00	0.00	0.00	0.50	0.50

比较转移矩阵$\mathbb{P}$和转移图，并确认它们包含的有关转移概率的信息一致。

为了找到从状态$i$转移到$j$的概率，只需找矩阵$i$行$j$列的值即可。

如果您知道起始状态，则可以使用$\mathbb{P}$找到任何有限路径的概率。例如，假设从1开始，那么它具有路径[2,2,3,4,3]的概率是

$$P(1, 2)P(2, 2)P(2, 3)P(3, 4)P(4, 3) \approx 0.4\%$$

0.5 * 0.5 * 0.25 * 0.25 * 0.25

0.00390625

MarkovChain对象的prob_of_path方法可以省去写乘法的麻烦。它将起始状态和路径的其余部分（在列表或数组中）作为其参数，并返回路径的概率。

reflecting_walk.prob_of_path(1, [2, 2, 3, 4, 3])

0.00390625

reflecting_walk.prob_of_path(1, [2, 2, 3, 4, 3, 5])

0.0

您可以使用simulate_path方法模拟链的路径。它有两个参数：起始状态和路径的步数。默认情况下，它返回一个由路径中的状态序列组成的数组。可选参数plot_path=True绘制模拟路径。运行几次下面的单元格，看看输出如何变化。

reflecting_walk.simulate_path(1, 7)

array([1, 2, 1, 2, 2, 2, 3, 2])

reflecting_walk.simulate_path(1, 10, plot_path=True)

$n$-阶转移矩阵

对于状态$i$和$j$, 耗费$n$步，从状态$i$转移为状态$j$的可能性，称为从$i$到$j$的$n$-阶转移概率。 形式上定义为：

$$P_n(i, j) ~ = ~ P(X_n = j \mid X_0 = i)$$

在这种表示方法中，一阶转移概率$P(i, j)$也可以写作$P_1(i, j)$。

$n$-阶转移概率$P_n(i, j)$可以使用$n$-阶转移矩阵$(i, j)$位置处的元素表示。对于任意状态$i$，$n$-阶转移矩阵的第$i$行包含了从状态$i$开始的链的条件分布$X_n$

MarkovChain的transition_matrix方法，使用$n$作为入参，并返回一个$n$-阶转移矩阵。以下是本节前面定义的粘性反转随机漫步的两阶转移矩阵。

reflecting_walk.transition_matrix(2)

	1	2	3	4	5
1	0.3750	0.5000	0.125	0.0000	0.0000
2	0.2500	0.4375	0.250	0.0625	0.0000
3	0.0625	0.2500	0.375	0.2500	0.0625
4	0.0000	0.0625	0.250	0.4375	0.2500
5	0.0000	0.0000	0.125	0.5000	0.3750

你可以轻松地手动计算各个条目。例如，$(1, 1)$是分两步从状态1进入状态1的可能性。有两种方法可以实现这一目标：

• [1, 1, 1]
• [1, 2, 1]

假设1是起始状态，则两条路径的总概率为$(0.5 \times 0.5) + (0.5 \times 0.25) = 0.375$.

由于马尔科夫性质，基于一阶转移概率，就能得到二阶转移概率。

一般而言，我们可以通过调节链条在时刻1时的位置，来计算$P_2(i, j)$

$$\begin{align*} P_2(i, j) ~ &= ~ P(X_2 = j \mid X_0 = i) \\ &= ~ \sum_k P(X_1 = k, X_2 = j \mid X_0 = i) \\ &= ~ \sum_k P(X_1 = k \mid X_0 = i)P(X_2 = j \mid X_1 = k) \\ &= ~ \sum_k P(i, k)P(k, j) \end{align*}$$

如上结果为$\mathbb{P} \times \mathbb{P} = \mathbb{P}^2$矩阵的$(i, j)$位置处元素。因此，二阶转移矩阵为$\mathbb{P}^2$。

通过归纳证明，能总结出，$n$-阶转移矩阵为$\mathbb{P}^n$。

这是粘性反转随机漫步的5步转移矩阵。

reflecting_walk.transition_matrix(5)

	1	2	3	4	5
1	0.246094	0.410156	0.234375	0.089844	0.019531
2	0.205078	0.363281	0.250000	0.136719	0.044922
3	0.117188	0.250000	0.265625	0.250000	0.117188
4	0.044922	0.136719	0.250000	0.363281	0.205078
5	0.019531	0.089844	0.234375	0.410156	0.246094

这是一个表示方法，但要使用矩阵，我们必须以Python识别为矩阵的形式表示它。方法get_transition_matrix为我们做到了这一点。需要步数$n$作为入参，并以numpy矩阵的格式返回$n$-阶转移矩阵。

对于粘性反转随机漫步，我们将从提取$\mathbb{P}$开始P作为矩阵refl_walk_P。

refl_walk_P = reflecting_walk.get_transition_matrix(1)
refl_walk_P

array([[ 0.5 ,  0.5 ,  0.  ,  0.  ,  0.  ],
       [ 0.25,  0.5 ,  0.25,  0.  ,  0.  ],
       [ 0.  ,  0.25,  0.5 ,  0.25,  0.  ],
       [ 0.  ,  0.  ,  0.25,  0.5 ,  0.25],
       [ 0.  ,  0.  ,  0.  ,  0.5 ,  0.5 ]])

让我们检查前面显示的5-阶转移矩阵是否与$\mathbb{P}^5$相同。您可以使用np.linalg.matrix_power将计算矩阵的非负整数次幂。第一个参数是矩阵，第二个参数是幂。

np.linalg.matrix_power(refl_walk_P, 5)

array([[ 0.24609375,  0.41015625,  0.234375  ,  0.08984375,  0.01953125],
       [ 0.20507812,  0.36328125,  0.25      ,  0.13671875,  0.04492188],
       [ 0.1171875 ,  0.25      ,  0.265625  ,  0.25      ,  0.1171875 ],
       [ 0.04492188,  0.13671875,  0.25      ,  0.36328125,  0.20507812],
       [ 0.01953125,  0.08984375,  0.234375  ,  0.41015625,  0.24609375]])

这确实与transition_matrix显示的矩阵相同，但难以阅读。

当我们想要在计算中使用$\mathbb{P}$，我们将使用此矩阵表示。对于显示和阅读，transition_matrix 则更好。

长期运行

要理解马尔科夫链的长跑行为，令$n$变大，并检查对于开始状态的每个$X_n$值。这些都包含在$n$-阶转移矩阵$\mathbb{P}^n$中。

如下展示了随机漫步中，$n = 25, 50$, 和 $100$情况下的$\mathbb{P}^n$ 。

reflecting_walk.transition_matrix(25)

	1	2	3	4	5
1	0.129772	0.256749	0.25	0.243251	0.120228
2	0.128374	0.254772	0.25	0.245228	0.121626
3	0.125000	0.250000	0.25	0.250000	0.125000
4	0.121626	0.245228	0.25	0.254772	0.128374
5	0.120228	0.243251	0.25	0.256749	0.129772

reflecting_walk.transition_matrix(50)

	1	2	3	4	5
1	0.125091	0.250129	0.25	0.249871	0.124909
2	0.125064	0.250091	0.25	0.249909	0.124936
3	0.125000	0.250000	0.25	0.250000	0.125000
4	0.124936	0.249909	0.25	0.250091	0.125064
5	0.124909	0.249871	0.25	0.250129	0.125091

reflecting_walk.transition_matrix(100)

	1	2	3	4	5
1	0.125	0.25	0.25	0.25	0.125
2	0.125	0.25	0.25	0.25	0.125
3	0.125	0.25	0.25	0.25	0.125
4	0.125	0.25	0.25	0.25	0.125
5	0.125	0.25	0.25	0.25	0.125

$\mathbb{P}^{100}$中每一行都一样！这意味着，对于随机漫步而言，在时刻100时的条件分布不依赖于其起始状态。链忘了他的起始点

你可以增加$n$，并且会看到$n$-阶转移矩阵保持一致。说明链已经达到稳定

稳定性是许多马尔可夫链的显着特性，也是本章的主题。

析构链

设$S$为一个有限状态或者可数无穷状态组成的集合。任一由集合$S$来索引行、列的随机矩阵，都是状态空间$S$下的某一马尔科夫链的转移矩阵。马尔可夫链的转移行为与矩阵变化保持一致。设置术语以讨论其中一些行为很有帮助。

连通

如果链可以从状态$i$转移到状态$j$，称之为$i$可达$j$，记作$i \rightarrow j$。通常，你能通过检查链的转移图来判定$i$ 是否可达$j$。一个$i \rightarrow j$正式的定义为：

• 存在一条转移路径，从$i$开始，到$j$结束。
• 等价的, 存在 $n > 0$ 使得 $P_n(i, j) > 0$。

当$i \rightarrow j$ 并且 $j \rightarrow i$时，称$i$连通$j$ 记作$i \leftrightarrow j$。

如果链的所有状态彼此连通，则该链被称为不可约。

上一节的粘性反转随机漫步是不可约的，因为链条可能从每个状态相互之间都是连通的。

周期

在离散时间工作存在缺陷。其中一点就是周期性。让我们从随机漫步的例子开始，其中每个步骤是基于公平硬币的投掷。假设从状态0开始。然后定义只能在如下时刻返回状态0：正面和反面出现的数量完全相等，因此投掷的数量必须是偶数。我们说状态0 有周期2

当链从状态$i$开始，并且经过$d$的倍数次数后能回到状态$i$，则称状态$i$ 有周期 $d$。$d$是所有能使 $P_n(i, i) > 0$的 $n$ 的最大公约数。

在上述描述的随机漫步中，所有的状态有周期2。

周期会导致长期行为的描述出现问题。例如：当状态$i$有周期3,序列$P_n(i, i)$可能看起来像”0, 0, positive, 0, 0, positive, $\ldots$”，因此限制声明可能会变得复杂。

在本课程中，我们将研究链条的长期行为，其中所有状态都是非周期性的，即它们具有周期1.换句话说，链条上无环。

你如何检查所有状态是否具有周期性？如果链是不可约的，那所有状态必须具有相同的周期。这个的证据并不困难，但我们不会这样做。因为这意味着如果一个链是不可约的，你就要找出其中每一个状态的周期，然后保证所有他状态都必须有相同的周期。

有些状态是很容易识别为非周期性的。如果1-阶转移概率$P(i, i)$ 为正，那么状态$i$是非周期性的。因为链可以保持在状态$i$任意长时间，其返回的结果是不成环的。

样例：析构链

考虑具有转移矩阵的链

	a	b	c	d	e
a	0	1	0	0	0
b	1	0	0	0	0
c	0	1/3	1/3	1/3	0
d	0	0	0	1/3	2/3
e	0	0	0	4/5	1/5

• 状态$a$和$b$相互连通，并且不和其他状态可达。因此称为连通类。小矩阵

	a	b
a	0	1
b	1	0

本身就是一个转移矩阵。尽管描述的是一个无聊的链，在状态$a$ 和$b$之间循环。$a$和$b$都有周期2。

• 状态$d$和$e$构成连通类，并且是非周期性的。

	d	e
d	1/3	2/3
e	4/5	1/5

• 状态$c$与自己连通，一旦转移为状态$b$或$d$, 就无法再返回。

在本课程中，我们将只使用有限状态空间上的不可约的非周期马尔可夫链。我们所说的大部分内容也适用于周期链，以及具有可数无限状态空间的链。

长期运行行为

每个具有有限状态空间的不可约和非周期性的马尔可夫链在状态转移一段时间后会表现出惊人的规律性。以下收敛定理的证明超出了本课程的范围，但您已经通过计算看到了结果。对于某些类别无限多个状态的马氏链，所有结果都更为正确。

收敛性

设$X_0, X_1, \ldots$ 是有限状态空间$S$上的不可约，非周期性的马尔科夫链。那么对于所有的状态$i$和$j$有

$$P_n(i, j) \to \pi(j) ~~~ \text{as } n \to \infty$$

换言之，对于$S$中任意的$i$和$j$，从$i$到$j$的$n$-阶转移概率会逼近一个极限，并且不依赖于$i$。
此外

• 对所有的状态$j$都有$\pi(j) > 0$ ，和

• $\sum_{j \in S} \pi(j) = 1$

也就是说，当$n \to \infty$， $n$-阶转移矩阵$\mathbb{P}^n$的中每一行的值都会等于同一个向量$\pi$，其中每一项都为正值。

限制特性

(i) 向量$\pi$是平衡方程 $\pi \mathbb{P} = \pi$的唯一解。

(ii) 如果对某些$n$，$X_n$的分布为$\pi$，那么，对于$m > n$，其分布 $X_m$ 也同样是$\pi$。因此，称$\pi$为链的静态或稳态分布。

(iii) 对于每一个状态$j$, 向量$\pi$的第$j$th 元素$\pi(j)$是链的长期值在$j$的预期。

我们假设收敛定理是正确的; 然后其他相关的特性推断起来就相对容易了。在本节的其余部分，我们将建立这些特性并查看它们的使用方式。

平衡方程

另$n \ge 0$，$i$和$j$是两个状态。然后

$$P_{n+1}(i, j) = \sum_{k \in S} P_n(i, k)P(k, j)$$

因此

$$\begin{align*} \lim_{n \to \infty} P_{n+1}(i, j) &= \lim_{n \to \infty} \sum_{k \in S} P_n(i, k)P(k, j) \\ \\ &= \sum_{k \in S} \big{(} \lim_{n \to \infty} P_n(i, k) \big{)} P(k, j) \end{align*}$$

因为$S$是有限的，我们可以交换极限与和。现在将收敛定理应用于平稳性：

$$\pi(j) = \sum_{k \in S} \pi(k)P(k, j)$$

这被称为平衡方程。

在矩阵表示方法中，如果你把$\pi$当做行向量，方程可以写为

$$\pi = \pi \mathbb{P} ~~~~~ \text{或者，等价的} ~~~~~ \pi\mathbb{P} = \pi$$

这有助于计算$\pi$时没有限制。

注意: 稳态不是状态空间$S$的元素。这是链条运行很长一段时间后的状况。让我们进一步研究这个问题。

平衡态和稳态

要想看看这些方程中的“平衡”是什么，就需要想象一下这个链的大量独立复制。例如，根据粘性反转随机漫步的转移概率，想象大量的粒子在状态1到5之间移动，并假设所有粒子在时刻1,2,3，……… 都彼此独立。

然后在任何时刻和任何状态$j$，有一些比例的粒子离开$j$，和另一些比例的粒子进入$j$。平衡方程表明这两个比例是相等的。

让我们通过再次查看方程来检查：对于任何状态$j$,

$$\pi(j) = \sum_{k \in S} \pi(k)P(k, j)$$

对于每一个$k \in S$ (包括 $k=j$)，以$\pi(k)$ 作为链运行很长一段时间后离开状态$k$的粒子的比例。等式左边是离开状态$j$的粒子的比例。等式右边求和的每一项都是离开状态$k$并转向状态$j$的粒子比例。求和之后，就是所有进入状态$j$的粒子。等式成立时，链是平衡的。

收敛于平稳性的定理表明，当$n$变大时，链趋于平衡。当链确实达到平衡时，对于$n$，其分布$X_n$为$\pi$，然后它保持平衡。原因：

$$P(X_{n+1} = j) = \sum_{i \in S} P(X_n = i)P(i, j) = \sum_{i \in S} \pi(i)P(i, j) = \pi(j)$$

通过平衡方程。现在使用归纳法。

特别是，如果链以其静止分布$\pi$开始，那么之后每一个$n$的分布$X_n$都是$\pi$。

唯一性

不难表明,如果平衡方程有解，那么它必须是$\pi$，$X_n$的边际分布的极限。我们不会做证明; 它基本上重复了我们用来推导平衡方程的步骤。你应该意识到，一个不可约的，非周期的，有限状态马尔可夫链只有一个稳态分布。

如果您碰巧猜测到平衡方程的解，这将特别有用。如果您猜到了一个概率分布解，那么您已经找到了链的稳态分布。

长期运行时各个状态占比

有状态$j$ ，令$I_m(j)$代表事件${X_m = j}$。链花费在状态$j$处转移次数比例，当转移次数从1至$n$，表述如下：

$$\frac{1}{n} \sum_{m=1}^n I_m(j)$$

因此，当链从状态$i$开始，链花费在状态$j$处转移次数比例预期为

$$\frac{1}{n} \sum_{m=1}^n E(I_m(j) \mid X_0 = i) = \frac{1}{n} \sum_{m=1}^n P(X_m = j \mid X_0 = i) = \frac{1}{n} \sum_{m=1}^n P_m(i, j)$$

现在回想一下实数序列的收敛性质：当$n \to \infty$时，$x_n \to x$，那么序列的均值也会收敛于$x$

$$\frac{1}{n} \sum_{m=1}^n x_m \to x ~~~ \text{as } n \to \infty$$

令$x_n = P_n(i, j)$。通过收敛性可得

$$P_n(i, j) \to \pi(j) ~~~ \text{as } n \to \infty$$

因此平均值也会收敛:

$$\frac{1}{n} \sum_{m=1}^n P_m(i, j) \to \pi(j) ~~~ \text{as } n \to \infty$$

因此，长期过程的链花费在状态$j$处转移次数比例预期为$\pi(j)$，其中$\pi$链的固定分布。

粘性反转漫步的稳态分布

我们在前面的部分对此进行了研究。转移图是

这是转移矩阵$\mathbb{P}$。

reflecting_walk

	1	2	3	4	5
1	0.50	0.50	0.00	0.00	0.00
2	0.25	0.50	0.25	0.00	0.00
3	0.00	0.25	0.50	0.25	0.00
4	0.00	0.00	0.25	0.50	0.25
5	0.00	0.00	0.00	0.50	0.50

MarkovChain的方法steady_state返回一个稳态分布$\pi$。 之前看到过这是$\mathbb{P}$行的极限。

reflecting_walk.steady_state()

Value	Probability
1	0.125
2	0.25
3	0.25
4	0.25
5	0.125

我们也可是使用平衡方程求解$\pi$。当然，这看起来有些多余，因为Python已经给出了$\pi$。当转移矩阵很大，且是分数的情况，使用Python是不错的操作。

根据平衡方程：

$$\pi(1) = \sum_{k=1}^s \pi(k)P(k, 1)$$

也就是说，使用$\pi$乘以$\mathbb{P}$中的1列

$$\pi(1) = \pi(1)\cdot 0.5 ~ + ~ \pi(2) \cdot 0.25 = 0.5\pi(1) + 0.25\pi(2)$$

按照相同的过程获得所有五个平衡方程：

$$\begin{align*} \pi(1) &= 0.5\pi(1) + 0.25\pi(2) \\ \pi(2) &= 0.5\pi(1) + 0.5\pi(2) + 0.25\pi(3) \\ \pi(3) &= 0.25\pi(2) + 0.5\pi(3) + 0.25\pi(4) \\ \pi(4) &= 0.25\pi(3) + 0.5\pi(4) + 0.5\pi(5) \\ \pi(5) &= 0.25\pi(4) + 0.5\pi(5) \end{align*}$$

一些观察结果使系统易于解决。

• 通过重新排列第一个等式，我们得到 $\pi(2) = 2\pi(1)$。
• 通过对称性, $\pi(1) = \pi(5)$ 和 $\pi(2) = \pi (4)$。
• 因为 $\pi(2) = \pi(4)$, 等式$\pi(3)$ 表明 $\pi(3) = \pi(2) = \pi(4)$。

所以$\pi$的分布为

$$\big{(} \pi(1), 2\pi(1), 2\pi(1), 2\pi(1), \pi(1) \big{)}$$

因为$\pi$是一个条件分布概率，其和为1。即 $8\pi(1)$的值为1，可以得到

$$\pi = \big{(} \frac{1}{8}, \frac{2}{8}, \frac{2}{8}, \frac{2}{8}, \frac{1}{8} \big{)}$$

这和我们用distribution计算$n=100$的结果一致。事实上，我们可以使用该方法steady_state来获得$\pi$:

这意味着从长远来看，这一部分的随机漫步预计将花费大约12.5％的时间在状态1，25％的时间用于状态2,3和4，其余12.5%的时间在状态5。

懒惰的随机循环漫步

现在让状态空间在圆上排列五个点。假设该过程从点1开始，并且在每个步骤中保持在概率为0.5的位置（因此是粘性的），或者移动到两个相邻点中的一个，每个概率为0.25，而不管其他移动。

换言之，除了$1 \rightarrow 5$ 和 $5 \rightarrow 1$这个漫步的转移和上面的随机漫步是相同的。 可以在转移图中总结此转移行为，请注意，所有状态的转移行为都是相同的。

在每一步中，下一步的动作都是通过从三个选项中随机选择和链的当前位置来确定的，而不是从它到达该位置的方式。所以这个过程就是马尔可夫链。我们称之为 $X_0, X_1, X_2, \ldots $ 并定义其转移矩阵。

s = np.arange(1, 6)

def circle_walk_probs(i, j):
        if i-j == 0:
            return 0.5
        elif abs(i-j) == 1:
            return 0.25
        elif abs(i-j) == 4:
            return 0.25
        else:
            return 0   
        
circle_walk = MarkovChain.from_transition_function(s, circle_walk_probs)

circle_walk

	1	2	3	4	5
1	0.50	0.25	0.00	0.00	0.25
2	0.25	0.50	0.25	0.00	0.00
3	0.00	0.25	0.50	0.25	0.00
4	0.00	0.00	0.25	0.50	0.25
5	0.25	0.00	0.00	0.25	0.50

由于转移行为的对称性，任何状态出现的概率，都不应该大于任何其他状态，因此$\pi(j)$的值是相同的。这可以用steady_state验证。

circle_walk.steady_state()

Value	Probability
1	0.2
2	0.2
3	0.2
4	0.2
5	0.2

样例

这里有两个例子来说明如何找到稳态分布以及如何使用它。

Ehrenfest的扩散模型

Paul Ehrenfest 提出了许多气体粒子扩散模型，其中一个我们将在这里研究。

该模型说有两个容器总共含有$N$个粒子。在每个瞬间，随机选择容器，并且独立于容器随机选择颗粒。然后将所选粒子放入所选容器中; 如果它已经在那个容器中，那就留在那里。

令$X_n$表示时刻$n$时容器1中的粒子数。那么$X_0, X_1, \ldots$是一个马尔科夫链，其转移概率描述如下：

\begin{equation}
P(i, j) =
\begin{cases}
\frac{N-i}{2N} & \text{if } j = i+1 \
\frac{1}{2} & \text{if } j = i \
\frac{i}{2N} & \text{if } j = i-1 \
0 & \text{otherwise}
\end{cases}
\end{equation}

这条链显然是不可约的。它是非周期性的，因为 $P(i, i) > 0$.

问题. 链的稳态分布是什么？

回答. 使用电脑， 所以，先找到$N=100$时的稳态分布，然后检查对于一般的$N$是否一致。

N = 100

states = np.arange(N+1)

def transition_probs(i, j):
    if j == i:
        return 1/2
    elif j == i+1:
        return (N-i)/(2*N)
    elif j == i-1:
        return i/(2*N)
    else:
        return 0

ehrenfest = MarkovChain.from_transition_function(states, transition_probs)
Plot(ehrenfest.steady_state(), edges=True)

这看起来很像二项式（100,1 / 2）分布。实际上它就是二项式（100,1 / 2）分布。既然你已经猜到了，你所要做的就是将它插入到平衡方程中并检查它们是否有效。

平衡方程是：

$$\begin{align*} \pi(0) &= \frac{1}{2}\pi(0) + \frac{1}{2N}\pi(1) \\ \pi(j) &= \frac{N-(j-1)}{2N}\pi(j-1) + \frac{1}{2}\pi(j) + \frac{j+1}{2N}\pi(j+1), ~~~ 1 \le j \le N-1 \\ \pi(N) &= \frac{1}{2N}\pi(N-1) + \frac{1}{2}\pi(N) \end{align*}$$

您已经通过查看$N=100$的结果猜测了答案。但是如果你想从头开始，你必须简化平衡方程并尝试用$\pi(0)$表示$\pi$的所有元素。你会得到：

$$\begin{align*} \pi(1) &= N\pi(0) \\ \\ \pi(2) &= \frac{N(N-1)}{2} \pi0 = \binom{N}{2} \pi(0) \end{align*}$$

可以归纳推导如下：

$$\pi(j) = \binom{N}{j} \pi(0)$$

换句话说，稳态分布分布与二项式系数成比例。所以当 $\pi(0) = 1/2^N$可以使所有元素的和为1,分布为二项分布 $(N, 1/2)$。

预期奖励

假设我长时间运行上一节中的懒惰反转随机漫步。如下，这是它的稳态分布。

stationary = reflecting_walk.steady_state()
stationary

Value	Probability
1	0.125
2	0.25
3	0.25
4	0.25
5	0.125

问题 1. 假设每次链条处于状态4时, 我赢得$$4$; 每次进入状态5, 我赢得$$5$; 否则我赢不到钱. 我奖励的期望是多少？

回答 1. 从长远来看，链条处于稳定状态。所以，有62.5%的概率，我赢不到钱，有25%的概率我赢$$4$，12.5%的概率，我赢$$5$。综上，计算可得奖励的期望为$$1.625$。

0*0.625 + 4*0.25 + 5*.125

1.625

问题 2. 假设每次链条处于状态$i$, 我抛$i$枚硬币，并记录正面次数。从长期来看，我每次得到正面个数的期望是多少？

回答 2. 每次链条处于状态$i$,期望得到$i/2$个正面。 当链处于稳态, 投掷硬币的期望数是3。所以，从长期来看，正面个数的期望是1.5。

stationary.ev()/2

1.5

这看上去是人为的，请考虑一下：假设我在上面玩游戏，并且在每一个动作中我告诉你我得到的头数，但我不告诉你链在哪个状态。我隐藏了潜在的马尔可夫链。如果您尝试重新创建马尔可夫链所采用的步骤序列，那么您正在使用隐马尔可夫模型。它们广泛用于模式识别，生物信息学和其他领域。