机器学习笔记之马尔可夫链蒙特卡洛方法——马尔可夫链与平稳分布
- 引言
- 回顾:蒙特卡洛方法
- 马尔可夫链与平稳分布
- 马尔可夫链
- 平稳分布
- 细致平衡
- 关于平稳分布的补充
- 马尔可夫链的本质
- 平稳分布的收敛性证明
- 相关总结
引言
上一节介绍了蒙特卡洛方法的具体思想 以及一些具体的蒙特卡洛方法。本节将介绍MCMC的另一个关键点:马尔可夫链(Markov Chain)。
回顾:蒙特卡洛方法
蒙特卡洛方法(Monte Carlo Method)的本质是通过采样的方式求解概率分布的期望信息。但主要问题在于采样过程中,由于待采样的概率分布 P ( X ) \mathcal P(\mathcal X) P(X)过于复杂(例如样本维度特征过高),从而无法实现 直接采样。
基于上述问题,介绍了两种采样方式:拒绝采样(Rejection Sampling)和重要性采样(Importance Sampling)。两种采样方式的共同点在于:
-
它们都没有直接对概率分布 P ( X ) \mathcal P(\mathcal X) P(X)进行采样,而是选择一个较容易采样的概率分布 Q ( X ) \mathcal Q(\mathcal X) Q(X)作为候选分布(Proposal Distribution),构建 P ( X ) \mathcal P(\mathcal X) P(X)和 Q ( X ) \mathcal Q(\mathcal X) Q(X)之间的关联关系,从而 通过采样 Q ( X ) \mathcal Q(\mathcal X) Q(X)来近似估计 P ( X ) \mathcal P(\mathcal X) P(X)的采样:
- 拒绝采样将某样本点 x ( i ) ∈ X x^{(i)} \in \mathcal X x(i)∈X在 P ( x ( i ) ) \mathcal P(x^{(i)}) P(x(i))和 M ⋅ Q ( x ( i ) ) \mathcal M \cdot \mathcal Q(x^{(i)}) M⋅Q(x(i))之间的差距关系映射成概率,从而基于该概率进行采样:
u ∼ U ( 0 , 1 ) u ⇔ ? P ( x ( i ) ) M ⋅ Q ( x ( i ) ) u \sim \mathcal U(0,1) \quad u \overset{\text{?}}{\Leftrightarrow} \frac{\mathcal P(x^{(i)})}{\mathcal M \cdot \mathcal Q(x^{(i)})} u∼U(0,1)u⇔?M⋅Q(x(i))P(x(i)) - 重要性采样引入候选分布 Q ( X ) \mathcal Q(\mathcal X) Q(X)将原始概率分布 P ( X ) \mathcal P(\mathcal X) P(X)的期望结果映射为带重要性权重的期望结果:
E P ( X ) [ f ( X ) ] = ∫ X [ P ( X ) Q ( X ) ⋅ f ( X ) ] d X = 1 N ∑ i = 1 N [ f ( x ( i ) ) ⋅ P ( x ( i ) ) Q ( x ( i ) ) ] \begin{aligned} \mathbb E_{\mathcal P(\mathcal X)} [f(\mathcal X)] & = \int_{\mathcal X} \left[\frac{\mathcal P(\mathcal X)}{\mathcal Q(\mathcal X)} \cdot f(\mathcal X)\right] d\mathcal X \\ & = \frac{1}{N} \sum_{i=1}^N \left[f(x^{(i)}) \cdot \frac{\mathcal P(x^{(i)})}{\mathcal Q(x^{(i)})}\right] \end{aligned} EP(X)[f(X)]=∫X[Q(X)P(X)⋅f(X)]dX=N1i=1∑N[f(x(i))⋅Q(x(i))P(x(i))]
- 拒绝采样将某样本点 x ( i ) ∈ X x^{(i)} \in \mathcal X x(i)∈X在 P ( x ( i ) ) \mathcal P(x^{(i)}) P(x(i))和 M ⋅ Q ( x ( i ) ) \mathcal M \cdot \mathcal Q(x^{(i)}) M⋅Q(x(i))之间的差距关系映射成概率,从而基于该概率进行采样:
-
在使用相应方法的过程中,对相应的候选分布 Q ( X ) \mathcal Q(\mathcal X) Q(X)存在限制条件:它们都需要候选分布 Q ( X ) \mathcal Q(\mathcal X) Q(X)与原始分布 P ( X ) \mathcal P(\mathcal X) P(X)相似:
- 拒绝采样需要满足 X \mathcal X X定义域内, M ⋅ Q ( X ) \mathcal M \cdot \mathcal Q(\mathcal X) M⋅Q(X)结果处处大于 P ( X ) \mathcal P(\mathcal X) P(X)的条件下,使 M \mathcal M M越小越好。但这不仅仅和 M \mathcal M M相关,选择的候选分布同样重要;
- 重要性采样需要选择的候选分布 Q ( X ) \mathcal Q(\mathcal X) Q(X)使重要性权重 P ( X ) Q ( X ) \frac{\mathcal P(\mathcal X)}{\mathcal Q(\mathcal X)} Q(X)P(X)尽可能地趋近于1。
P ( X ) Q ( X ) ⇔ 1 \frac{\mathcal P(\mathcal X)}{\mathcal Q(\mathcal X)} \Leftrightarrow 1 Q(X)P(X)⇔1
实际情况是,在基于高维空间中的概率分布 P ( X ) \mathcal P(\mathcal X) P(X),很难找到一个与 P ( X ) \mathcal P(\mathcal X) P(X)近似并且采样简单的概率分布 Q ( X ) \mathcal Q(\mathcal X) Q(X)。因此,需要观察是否存在其他的采样方法——马尔可夫链蒙特卡洛方法(Markov Chain Monte Carlo,MCMC)。
马尔可夫链与平稳分布
马尔可夫链
马尔可夫链是一种具有马尔可夫性质,并且存在于离散的时间与状态下的随机过程。
随机过程研究的对象不同于以往的随机变量,而是随机变量序列。
- 通常情况下,随机变量 X \mathcal X X可能服从某一概率分布 P ( X ) P(\mathcal X) P(X):
X ∼ P ( X ) \mathcal X \sim P(\mathcal X) X∼P(X) - 随机过程研究的对象可能是基于时间或者空间变化的若干离散随机变量组成的序列。通常用 { X T } \{\mathcal X_T\} {XT}表示:
有可能是有限的,也有可能是无限的~
{ X T } = { X 1 , X 2 , ⋯ , X t , X t + 1 , ⋯ } \{\mathcal X_T\} = \{\mathcal X_1,\mathcal X_2,\cdots,\mathcal X_t,\mathcal X_{t+1},\cdots\} {XT}={X1,X2,⋯,Xt,Xt+1,⋯}
而马尔可夫链具备马尔可夫性质。即隐马尔可夫模型的背景介绍中提到的齐次马尔可夫假设:给定当前状态随机变量 X t \mathcal X_t Xt的条件下,未来状态的随机变量与过去状态无关。以一阶齐次马尔可夫假设为例,数学符号表示如下:
P ( X t + 1 ∣ X 1 , X 2 , ⋯ , X t ) = P ( X t + 1 ∣ X t ) P(\mathcal X_{t+1} \mid \mathcal X_1,\mathcal X_2,\cdots,\mathcal X_{t}) = P(\mathcal X_{t+1} \mid \mathcal X_t) P(Xt+1∣X1,X2,⋯,Xt)=P(Xt+1∣Xt)
由于各状态的随机变量内部均是离散的。因此,和隐马尔可夫模型类似,定义 A \mathcal A A为状态转移矩阵。
- 基于离散的随机变量,定义任一状态的随机变量 X \mathcal X X存在 K \mathcal K K种取法。即:
X 1 , ⋯ , X T ∈ { x 1 , x 2 , ⋯ , x K } \mathcal X_1,\cdots,\mathcal X_T \in \{x_1,x_2,\cdots,x_{\mathcal K}\} X1,⋯,XT∈{x1,x2,⋯,xK} - 因而,状态转移矩阵是一个 K × K \mathcal K\times \mathcal K K×K的方阵:
A = [ a i j ] K × K ( i , j ∈ { 1 , 2 , ⋯ , K } ) \mathcal A = [a_{ij}]_{\mathcal K\times \mathcal K} \quad (i,j \in \{1,2,\cdots,\mathcal K\}) A=[aij]K×K(i,j∈{1,2,⋯,K})
方阵 A \mathcal A A中的每一个元素 a i j a_{ij} aij均可表示如下:
一个关于‘马尔可夫性质’的条件概率。
a i j = P ( X t + 1 = x j ∣ X t = x i ) a_{ij} = P(\mathcal X_{t+1}=x_j \mid \mathcal X_{t} = x_i) aij=P(Xt+1=xj∣Xt=xi)
平稳分布
基于上述介绍,一个马尔可夫链的概率图 表示如下:
针对概率图中的任一随机变量 X i ( i = 1 , 2 ⋯ ) \mathcal X_i (i=1,2\cdots) Xi(i=1,2⋯),它的选择结果是离散的,因而每一个随机变量均对应一个概率分布。这里以 X 1 \mathcal X_1 X1为例,构建它的 概率分布:
- 随机变量 X 1 \mathcal X_1 X1可选择的离散结果以及对应概率表示如下:
X 1 \mathcal X_1 X1可选择结果 | x 1 x_1 x1 | x 2 x_2 x2 | ⋯ \cdots ⋯ | x K x_{\mathcal K} xK |
---|---|---|---|---|
各结果对应概率 | Q ( X 1 = x 1 ) \mathcal Q(\mathcal X_1 = x_1) Q(X1=x1) | Q ( X 1 = x 2 ) \mathcal Q(\mathcal X_1 = x_2) Q(X1=x2) | ⋯ \cdots ⋯ | Q ( X 1 = x K ) \mathcal Q(\mathcal X_1 = x_{\mathcal K}) Q(X1=xK) |
- 根据概率密度积分公式,自然有:
∑ i = 1 K Q ( X 1 = x i ) = 1 \sum_{i=1}^{\mathcal K}\mathcal Q(\mathcal X_1 = x_i) = 1 i=1∑KQ(X1=xi)=1 - 因此,定义 π ( X 1 ) \pi(\mathcal X_1) π(X1)为随机变量 X 1 \mathcal X_1 X1的概率分布:
π ( X 1 ) = { Q ( X 1 = x 1 ) , Q ( X 1 = x 2 ) ⋯ , Q ( X 1 = x K ) } \pi(\mathcal X_1) = \{\mathcal Q(\mathcal X_1 = x_1),\mathcal Q(\mathcal X_1 = x_2)\cdots,\mathcal Q(\mathcal X_1 = x_{\mathcal K})\} π(X1)={Q(X1=x1),Q(X1=x2)⋯,Q(X1=xK)} - 同理,其他随机变量同样存在 相应的概率分布:
π ( X 1 ) , π ( X 2 ) , ⋯ π ( X t ) , π ( X t + 1 ) , ⋯ \pi(\mathcal X_1),\pi(\mathcal X_2),\cdots\pi(\mathcal X_t),\pi(\mathcal X_{t+1}),\cdots π(X1),π(X2),⋯π(Xt),π(Xt+1),⋯
而平稳分布(Stationary Distribution)表示随着时间状态的推移,基于马尔可夫性质,最终会收敛至一个相对稳定的分布。
假设状态转移过程中到达第 m m m时刻 实现平稳分布,即 m m m时刻之后的所有概率分布均可视为相同的概率分布。即:
π ( X m ) = π ( X m + 1 ) = ⋯ \pi(\mathcal X_m) = \pi(\mathcal X_{m+1}) = \cdots π(Xm)=π(Xm+1)=⋯
以随机变量 X t \mathcal X_t Xt向随机变量 X t + 1 \mathcal X_{t+1} Xt+1的状态转移过程 为例,基于马尔可夫性质,概率分布之间的关联关系表示如下:
定义一个符号:
π ( X t + 1 = x ∗ ) \pi(\mathcal X_{t+1} = x^*) π(Xt+1=x∗):表示在t+1时刻的概率分布
π ( X t + 1 ) \pi(\mathcal X_{t+1}) π(Xt+1)下,随机变量
X t + 1 \mathcal X_{t+1} Xt+1选择一个具体结果
x ∗ x^* x∗;
π ( X t + 1 = x ∗ ) = ∑ x 1 , ⋯ , x K π ( X t = x ) ⋅ P ( X t + 1 = x ∗ ∣ X t = x ) x ∈ { x 1 , x 2 , ⋯ , x K } \pi(\mathcal X_{t+1} = x^*) = \sum_{x_1,\cdots,x_{\mathcal K}} \pi(\mathcal X_t = x) \cdot P(\mathcal X_{t+1} = x^* \mid \mathcal X_t = x) \quad x \in \{x_1,x_2,\cdots,x_{\mathcal K}\} π(Xt+1=x∗)=x1,⋯,xK∑π(Xt=x)⋅P(Xt+1=x∗∣Xt=x)x∈{x1,x2,⋯,xK}
假设存在一个概率分布 π ∗ ( X ) \pi^*(\mathcal X) π∗(X):
π ∗ ( X ) = { π ∗ ( X = x 1 ) , π ∗ ( X = x 2 ) , ⋯ , π ∗ ( X = x K ) } ∑ i = 1 K π ∗ ( X = x i ) = 1 \pi^*(\mathcal X) = \{\pi^*(\mathcal X = x_1),\pi^*(\mathcal X = x_2),\cdots,\pi^*(\mathcal X = x_{\mathcal K})\} \\ \sum_{i=1}^{\mathcal K} \pi^*(\mathcal X = x_i) = 1 π∗(X)={π∗(X=x1),π∗(X=x2),⋯,π∗(X=xK)}i=1∑Kπ∗(X=xi)=1
该概率分布在任意时刻都能满足:
π ∗ ( X = x ∗ ) = ∑ x 1 , ⋯ , x K π ∗ ( X = x ) ⋅ P ( X = x ∗ ∣ X = x ) \pi^*(\mathcal X = x^*) = \sum_{x_1,\cdots,x_{\mathcal K}} \pi^*(\mathcal X = x) \cdot P(\mathcal X = x^* \mid \mathcal X = x) π∗(X=x∗)=x1,⋯,xK∑π∗(X=x)⋅P(X=x∗∣X=x)
则 π ∗ ( X ) \pi^*(\mathcal X) π∗(X)被称为平稳分布。通俗点说就是分布 π ∗ ( X ) \pi^*(\mathcal X) π∗(X)在任意时刻都能满足齐次马尔可夫假设。从而其他时刻随机变量的概率分布 都可以使用平稳分布进行表示:
π ( X 1 ) , π ( X 2 ) , ⋯ π ( X t ) , π ( X t + 1 ) , ⋯ → π ∗ ( X ) \pi(\mathcal X_1),\pi(\mathcal X_2),\cdots\pi(\mathcal X_t),\pi(\mathcal X_{t+1}),\cdots \to \pi^*(\mathcal X) π(X1),π(X2),⋯π(Xt),π(Xt+1),⋯→π∗(X)
细致平衡
平稳分布 π ∗ ( X ) \pi^*(\mathcal X) π∗(X)在 马尔可夫链的任意时刻均可以进行表示。如果我们将得到平稳分布 作为目标,去构建马尔科夫链。满足什么条件的概率分布才是平稳分布?
这里介绍一个条件:细致平衡(Detail Balance)。
细致平衡是一个统计物理学名词,所谓热力学平衡是通过细致平衡达到的。而在马尔可夫过程中,它表示一种可逆性条件。具体公式表示如下:
π ( X = x i ) ⋅ a i j = π ( X = x j ) ⋅ a j i i , j ∈ { 1 , 2 , ⋯ , K } \pi(\mathcal X = x_i) \cdot a_{ij} = \pi(\mathcal X = x_j) \cdot a_{ji} \quad i,j \in \{1,2,\cdots,\mathcal K\} π(X=xi)⋅aij=π(X=xj)⋅ajii,j∈{1,2,⋯,K}
对应的图形关系如下所示:
如果将细致平衡延伸至更泛化的条件,具体公式表示如下:
其中
x x x可以表示
{ x 1 , ⋯ , x K } \{x_1,\cdots,x_{\mathcal K}\} {x1,⋯,xK}中的任意一个元素,
x ∗ x^* x∗表示经过状态转移后产生的未来状态的某元素。
π ( X = x ) ⋅ P ( X = x ∗ ∣ X = x ) = π ( X = x ∗ ) ⋅ P ( X = x ∣ X = x ∗ ) \pi(\mathcal X = x) \cdot P(\mathcal X = x^* \mid \mathcal X = x) = \pi(\mathcal X = x^*) \cdot P(\mathcal X = x \mid \mathcal X = x^*) π(X=x)⋅P(X=x∗∣X=x)=π(X=x∗)⋅P(X=x∣X=x∗)
细致平衡是马尔科夫链存在平稳分布的充分非必要条件。即某一时刻如果马尔可夫链中的各元素满足细致平衡条件,那么该时刻随机变量 X \mathcal X X的概率分布就是平稳分布。
- 首先观察平稳分布的关系式:
x x x项既在
π ∗ ( X = x ) \pi^*(\mathcal X = x) π∗(X=x)又在
P ( X = x ∗ ∣ X = x ) P(\mathcal X = x^* \mid \mathcal X = x) P(X=x∗∣X=x)中,无法进行化简;
∑ x 1 , ⋯ , x K π ∗ ( X = x ) ⋅ P ( X = x ∗ ∣ X = x ) d x \sum_{x_1,\cdots,x_{\mathcal K}} \pi^*(\mathcal X = x) \cdot P(\mathcal X = x^* \mid \mathcal X = x) dx x1,⋯,xK∑π∗(X=x)⋅P(X=x∗∣X=x)dx - 将细致平衡条件带入,即使用 π ∗ ( X = x ∗ ) ⋅ P ( X = x ∣ X = x ∗ ) \pi^*(\mathcal X = x^*) \cdot P(\mathcal X = x \mid \mathcal X = x^*) π∗(X=x∗)⋅P(X=x∣X=x∗)替换上式中的 π ∗ ( X = x ) ⋅ P ( X = x ∗ ∣ X = x ) \pi^*(\mathcal X = x) \cdot P(\mathcal X = x^* \mid \mathcal X = x) π∗(X=x)⋅P(X=x∗∣X=x):
∑ x 1 , ⋯ , x K π ∗ ( X = x ∗ ) ⋅ P ( X = x ∣ X = x ∗ ) d x = π ∗ ( X = x ∗ ) ⋅ ∑ x 1 , ⋯ , x K P ( X = x ∣ X = x ∗ ) d x = π ∗ ( X = x ∗ ) \begin{aligned} & \sum_{x_1,\cdots,x_{\mathcal K}} \pi^*(\mathcal X = x^*) \cdot P(\mathcal X = x \mid \mathcal X = x^*) dx \\ & = \pi^*(\mathcal X = x^*) \cdot \sum_{x_1,\cdots,x_{\mathcal K}} P(\mathcal X = x \mid \mathcal X = x^*) dx \\ & = \pi^*(\mathcal X = x^*) \end{aligned} x1,⋯,xK∑π∗(X=x∗)⋅P(X=x∣X=x∗)dx=π∗(X=x∗)⋅x1,⋯,xK∑P(X=x∣X=x∗)dx=π∗(X=x∗)
至此,通过细致平衡条件,可以推出平稳分布的关系式。我们在构建马尔可夫链时,只需要满足细致平衡条件,就可以得到平稳分布,而不是通过定义的条件去求解平稳分布。
关于平稳分布的补充
马尔可夫链的本质
重新观察平稳分布的关系式:
π ∗ ( X = x ∗ ) = ∑ x 1 , ⋯ , x K π ∗ ( X = x ) ⋅ P ( X = x ∗ ∣ X = x ) \pi^*(\mathcal X = x^*) = \sum_{x_1,\cdots,x_{\mathcal K}} \pi^*(\mathcal X = x) \cdot P(\mathcal X = x^* \mid \mathcal X = x) π∗(X=x∗)=x1,⋯,xK∑π∗(X=x)⋅P(X=x∗∣X=x)
需要注意的点是:无论是 π ∗ ( X = x ) \pi^*(\mathcal X = x) π∗(X=x)还是 π ∗ ( X = x ∗ ) \pi^*(\mathcal X = x^*) π∗(X=x∗),它们都是计算过程中产生的中间结果,每个状态的转移过程是基于概率的随机选择,不受人为控制。
而真正固定不变的量是状态转移矩阵。通过构建有效的状态转移矩阵可以约束每个状态的元素向下一状态转移的过程。
最终,通过状态转移矩阵的约束,随着状态的变化,使每个状态的分布逼近至平稳分布。
平稳分布的收敛性证明
-
详细观察状态转移矩阵 A \mathcal A A,将其展开:
A = ( a 11 , a 12 , ⋯ , a 1 K a 21 , a 22 , ⋯ , a 2 K ⋮ a K 1 , a K 2 , ⋯ , a K K ) K × K a i j = P ( X t + 1 = x j ∣ X t = x i ) \mathcal A = \begin{pmatrix} a_{11},a_{12},\cdots,a_{1 \mathcal K} \\ a_{21},a_{22},\cdots,a_{2 \mathcal K} \\ \vdots \\ a_{\mathcal K 1},a_{\mathcal K 2},\cdots,a_{\mathcal K \mathcal K}\end{pmatrix}_{\mathcal K \times \mathcal K} \quad a_{ij} = P(\mathcal X_{t+1} = x_j \mid \mathcal X_{t} = x_i) A=⎝ ⎛a11,a12,⋯,a1Ka21,a22,⋯,a2K⋮aK1,aK2,⋯,aKK⎠ ⎞K×Kaij=P(Xt+1=xj∣Xt=xi)
首先,矩阵 A \mathcal A A中的每一行结果之和为 1 1 1。以第一行为例,即:
其表示的具体意义是:在当前状态选择
x 1 x_1 x1的条件下,必然要转移至下一状态。而下一状态只包含离散的
K \mathcal K K种,只能从
K \mathcal K K种转移结果中选择一个,无其他选择;
如果随机变量是连续的 -> 将描述修改为‘范围’即可。
∑ j = 1 K a 1 j = 1 \sum_{j=1}^{\mathcal K} a_{1j} = 1 j=1∑Ka1j=1
其次,矩阵 A \mathcal A A中的每一列结果之和为1。以第一列为例,即:
其表示的具体意义是:确定了状态转移的目标是
x 1 x_1 x1的条件下,同样也只能从离散的
K \mathcal K K种状态中进行选择,并进行转移。
∑ i = 1 K a i 1 = 1 \sum_{i=1}^{\mathcal K} a_{i1} = 1 i=1∑Kai1=1
我们也称 A \mathcal A A为随机矩阵(Stochastic Matrix) -
基于上述随机矩阵的逻辑,从 t t t时刻开始转移,转移至 t + 1 t+1 t+1时刻选择 x j x_j xj的概率如何表示?
实际上,并不知道
t t t时刻到底是哪一项元素
x i ( i = 1 , 2 , ⋯ , K ) x_i(i=1,2,\cdots,\mathcal K) xi(i=1,2,⋯,K),而该元素的产生是由
t t t时刻状态的概率分布决定的。
π ( X t + 1 = x j ) = ∑ i = 1 K π ( X t = x i ) ⋅ a i j \pi(\mathcal X_{t+1} = x_j) = \sum_{i=1}^{\mathcal K} \pi(\mathcal X_t = x_i) \cdot a_{ij} π(Xt+1=xj)=i=1∑Kπ(Xt=xi)⋅aij
观察上述公式,包含求和符号,并且包含随机矩阵 A \mathcal A A的元素项 a i j a_{ij} aij。尝试将上式用矩阵表达出来:- 定义 t + 1 t+1 t+1时刻状态下随机变量 X t + 1 \mathcal X_{t+1} Xt+1的概率分布表示为:
π ( X t + 1 ) = [ π ( X t + 1 = x 1 ) , π ( X t + 1 = x 2 ) , ⋯ , π ( X t + 1 = x K ) ] 1 × K \pi(\mathcal X_{t+1}) = \left[\pi(\mathcal X_{t+1} = x_1),\pi(\mathcal X_{t+1} = x_2),\cdots,\pi(\mathcal X_{t+1} = x_{\mathcal K})\right]_{1 \times \mathcal K} π(Xt+1)=[π(Xt+1=x1),π(Xt+1=x2),⋯,π(Xt+1=xK)]1×K - 基于上述定义,利用上式对 π ( X t + 1 ) \pi(\mathcal X_{t+1}) π(Xt+1)进行表示:
写成列向量的形式,便于观察~
π ( X t + 1 ) = [ ∑ i = 1 K π ( X t = x i ) ⋅ a i 1 ∑ i = 2 K π ( X t = x i ) ⋅ a i 2 ⋮ ∑ i = 1 K π ( X t = x i ) ⋅ a i K ] 1 × K T \pi(\mathcal X_{t+1}) = \begin{bmatrix} \sum_{i=1}^{\mathcal K}\pi(\mathcal X_t = x_i) \cdot a_{i1} \\ \sum_{i=2}^{\mathcal K}\pi(\mathcal X_t = x_i) \cdot a_{i2} \\ \vdots \\ \sum_{i=1}^{\mathcal K} \pi(\mathcal X_t = x_i) \cdot a_{i\mathcal K}\end{bmatrix}^{T}_{1 \times \mathcal K} π(Xt+1)=⎣ ⎡∑i=1Kπ(Xt=xi)⋅ai1∑i=2Kπ(Xt=xi)⋅ai2⋮∑i=1Kπ(Xt=xi)⋅aiK⎦ ⎤1×KT
以第一项 ∑ i = 1 K π ( X = x i ) ⋅ a i 1 \sum_{i=1}^{\mathcal K} \pi(\mathcal X = x_i) \cdot a_{i1} ∑i=1Kπ(X=xi)⋅ai1为例:
∑ i = 1 K π ( X = x i ) ⋅ a i 1 = π ( X t = x 1 ) ⋅ a 11 + π ( X t = x 2 ) ⋅ a 21 + ⋯ + π ( X t = x K ) ⋅ a K 1 = [ π ( X t = x 1 ) , ⋯ , π ( X t = x K ) ] ( a 11 a 21 ⋮ a K 1 ) \begin{aligned} \sum_{i=1}^{\mathcal K} \pi(\mathcal X = x_i) \cdot a_{i1} & = \pi(\mathcal X_t = x_1) \cdot a_{11} + \pi(\mathcal X_t = x_2) \cdot a_{21} +\cdots + \pi(\mathcal X_t = x_{\mathcal K}) \cdot a_{\mathcal K 1} \\ & = [\pi(\mathcal X_t = x_1),\cdots,\pi(\mathcal X_t = x_{\mathcal K})] \begin{pmatrix}a_{11} \\ a_{21} \\ \vdots \\ a_{\mathcal K 1}\end{pmatrix} \end{aligned} i=1∑Kπ(X=xi)⋅ai1=π(Xt=x1)⋅a11+π(Xt=x2)⋅a21+⋯+π(Xt=xK)⋅aK1=[π(Xt=x1),⋯,π(Xt=xK)]⎝ ⎛a11a21⋮aK1⎠ ⎞
同理,其他项也可以表示为上述形式:
[ π ( X t = x 1 ) , ⋯ , π ( X t = x K ) ] ( a 12 a 22 ⋮ a K 2 ) , ⋯ , [ π ( X t = x 1 ) , ⋯ , π ( X t = x K ) ] ( a 1 K a 2 K ⋮ a K K ) [\pi(\mathcal X_t = x_1),\cdots,\pi(\mathcal X_t = x_{\mathcal K})] \begin{pmatrix}a_{12} \\ a_{22} \\ \vdots \\ a_{\mathcal K 2}\end{pmatrix},\cdots, [\pi(\mathcal X_t = x_1),\cdots,\pi(\mathcal X_t = x_{\mathcal K})] \begin{pmatrix}a_{1\mathcal K} \\ a_{2\mathcal K} \\ \vdots \\ a_{\mathcal K \mathcal K}\end{pmatrix} [π(Xt=x1),⋯,π(Xt=xK)]⎝ ⎛a12a22⋮aK2⎠ ⎞,⋯,[π(Xt=x1),⋯,π(Xt=xK)]⎝ ⎛a1Ka2K⋮aKK⎠ ⎞
将上述所有项进行合并,有:
[ π ( X t = x 1 ) , ⋯ , π ( X t = x K ) ] [\pi(\mathcal X_t = x_1),\cdots,\pi(\mathcal X_t = x_{\mathcal K})] [π(Xt=x1),⋯,π(Xt=xK)]就是
t t t时刻状态下随机变量
X t \mathcal X_t Xt的概率分布
π ( X t ) \pi(\mathcal X_t) π(Xt);
π ( X t + 1 ) = [ π ( X t = x 1 ) , ⋯ , π ( X t = x K ) ] ⋅ ( a 11 , a 12 , ⋯ , a 1 K a 21 , a 22 , ⋯ , a 2 K ⋮ a K 1 , a K 2 , ⋯ , a K K ) K × K = π ( X t ) ⋅ A \begin{aligned} \pi(\mathcal X_{t+1}) & = [\pi(\mathcal X_t = x_1),\cdots,\pi(\mathcal X_t = x_{\mathcal K})] \cdot \begin{pmatrix} a_{11},a_{12},\cdots,a_{1\mathcal K} \\ a_{21},a_{22},\cdots,a_{2\mathcal K} \\ \vdots \\ a_{\mathcal K 1},a_{\mathcal K 2},\cdots, a_{\mathcal K \mathcal K} \end{pmatrix}_{\mathcal K \times \mathcal K} \\ & = \pi(\mathcal X_t) \cdot \mathcal A \end{aligned} π(Xt+1)=[π(Xt=x1),⋯,π(Xt=xK)]⋅⎝ ⎛a11,a12,⋯,a1Ka21,a22,⋯,a2K⋮aK1,aK2,⋯,aKK⎠ ⎞K×K=π(Xt)⋅A - 定义 t + 1 t+1 t+1时刻状态下随机变量 X t + 1 \mathcal X_{t+1} Xt+1的概率分布表示为:
-
至此,找到了 t t t时刻与 t + 1 t+1 t+1时刻随机变量概率分布 之间的关联关系。实际上,从初始概率分布 π ( X i n i t ) \pi(\mathcal X_{init}) π(Xinit)开始,每一时间节点均乘一次随机矩阵 A \mathcal A A:
π ( X t + 1 ) = π ( X i n i t ) ⋅ A t \pi(\mathcal X_{t+1}) = \pi(\mathcal X_{init}) \cdot \mathcal A^{t} π(Xt+1)=π(Xinit)⋅At
根据随机矩阵 A \mathcal A A的一条特殊性质:随机矩阵各特征值的绝对值 ≤ 1 \leq 1 ≤1。
基于该性质,对 A \mathcal A A进行特征值分解:
A = P Λ P − 1 \mathcal A = \mathcal P \Lambda\mathcal P^{-1} A=PΛP−1
而对角矩阵 Λ \Lambda Λ对角线上元素即 A \mathcal A A的特征值结果:
Λ = ( λ 1 , 0 , ⋯ , 0 0 , λ 2 , ⋯ , 0 ⋮ 0 , 0 , ⋯ , λ K ) λ i ≤ 1 ; i = 1 , 2 , ⋯ , K \Lambda = \begin{pmatrix}\lambda_1 ,0,\cdots, 0 \\ 0,\lambda_2,\cdots,0 \\ \vdots \\0,0, \cdots,\lambda_{\mathcal K}\end{pmatrix} \quad \lambda_i \leq1;i = 1,2,\cdots,\mathcal K Λ=⎝ ⎛λ1,0,⋯,00,λ2,⋯,0⋮0,0,⋯,λK⎠ ⎞λi≤1;i=1,2,⋯,K
将 A \mathcal A A的分解结果带入 π ( X t + 1 ) \pi(\mathcal X_{t+1}) π(Xt+1)中:
π ( X t + 1 ) = π ( X i n i t ) ⋅ A t = π ( X i n i t ) ( P Λ P − 1 ) t = π ( X i n i t ) ( P Λ P − 1 P Λ P − 1 ⋯ P Λ P − 1 ) \begin{aligned}\pi(\mathcal X_{t+1}) & = \pi(\mathcal X_{init}) \cdot \mathcal A^t \\ & = \pi(\mathcal X_{init}) (\mathcal P\Lambda\mathcal P^{-1})^t \\ & = \pi(\mathcal X_{init}) (\mathcal P\Lambda\mathcal P^{-1}\mathcal P\Lambda\mathcal P^{-1}\cdots\mathcal P\Lambda\mathcal P^{-1}) \end{aligned} π(Xt+1)=π(Xinit)⋅At=π(Xinit)(PΛP−1)t=π(Xinit)(PΛP−1PΛP−1⋯PΛP−1)
由于 P − 1 P = E \mathcal P^{-1}\mathcal P = E P−1P=E,即单位矩阵。因此, π ( X t + 1 ) \pi(\mathcal X_{t+1}) π(Xt+1)的结果表示如下:
π ( X t + 1 ) = π ( X i n i t ) ⋅ P Λ t P − 1 = π ( X i n i t ) ⋅ P ( λ 1 t , 0 , ⋯ , 0 0 , λ 2 t , ⋯ , 0 ⋮ 0 , 0 , ⋯ , λ K t ) P − 1 \begin{aligned} \pi(\mathcal X_{t+1}) & = \pi(\mathcal X_{init}) \cdot \mathcal P \Lambda^t \mathcal P^{-1} \\ & = \pi(\mathcal X_{init}) \cdot \mathcal P \begin{pmatrix}\lambda_1^t ,0,\cdots, 0 \\ 0,\lambda_2^t,\cdots,0 \\ \vdots \\0,0, \cdots,\lambda_{\mathcal K}^t\end{pmatrix}\mathcal P^{-1} \end{aligned} π(Xt+1)=π(Xinit)⋅PΛtP−1=π(Xinit)⋅P⎝ ⎛λ1t,0,⋯,00,λ2t,⋯,0⋮0,0,⋯,λKt⎠ ⎞P−1 -
当转移次数足够多的时候,即 t t t足够大的时候,由于 λ 1 , ⋯ , λ K \lambda_1,\cdots,\lambda_{\mathcal K} λ1,⋯,λK均小于等于1恒成立,则有:
lim t → ∞ λ i t = 0 ( i = 1 , 2 , ⋯ , K ) \mathop{\lim}\limits_{t \to \infty} \lambda_i^t = 0 \quad (i=1,2,\cdots,\mathcal K) t→∞limλit=0(i=1,2,⋯,K)
假设转移次数分别执行了 m m m次和 m + 1 m+1 m+1次,它们对应随机变量的概率分布表示如下:
π ( X m + 1 ) = π ( X i n i t ) ⋅ P Λ m P − 1 π ( X m + 2 ) = π ( X i n i t ) ⋅ P Λ ( m + 1 ) P − 1 \pi(\mathcal X_{m+1}) = \pi(\mathcal X_{init}) \cdot \mathcal P\Lambda^{m}\mathcal P^{-1} \\ \pi(\mathcal X_{m+2}) = \pi(\mathcal X_{init}) \cdot \mathcal P\Lambda^{(m+1)} \mathcal P^{-1} π(Xm+1)=π(Xinit)⋅PΛmP−1π(Xm+2)=π(Xinit)⋅PΛ(m+1)P−1
而此时的 对角阵 Λ m , λ m + 1 \Lambda^{m},\lambda^{m+1} Λm,λm+1中的各元素均无限趋近于0(特征值为1的结果除外),因此,在 m m m足够大时:
Λ m = Λ m + 1 \Lambda^{m} = \Lambda^{m+1} Λm=Λm+1
从而有:
π ( X m + 1 ) = π ( X m + 2 ) \pi(\mathcal X_{m+1}) = \pi(\mathcal X_{m+2}) π(Xm+1)=π(Xm+2) -
基于上述推导过程,如果转移次数 t > m t > m t>m时,有:
π ( X e n d ) \pi(\mathcal X_{end}) π(Xend)表示转移过程的终结状态~
π ( X m + 1 ) = π ( X m + 2 ) = ⋯ = π ( X e n d ) \pi(\mathcal X_{m+1}) = \pi(\mathcal X_{m+2}) = \cdots = \pi(\mathcal X_{end}) π(Xm+1)=π(Xm+2)=⋯=π(Xend)
相关总结
- 平稳分布是马尔可夫链的一种性质,这种性质是基于状态转移矩阵是 双随机矩阵 的性质得到,和各时刻随机变量的概率分布之间没有关系;
双随机矩阵即‘矩阵各行,列结果之和为1’的矩阵。
- 细致平衡是马尔可夫链存在平稳分布 的充分非必要条件,因此它是衡量平稳分布的重要指标;
- 对于马尔可夫链,唯一能够调整(按照人的意志)是状态转移矩阵,换句话说,随机变量 X \mathcal X X的状态转移的过程是基于概率的随机结果,不因人的意志变化而变化,我们只能通过调整状态转移概率来影响模型。
下一节将介绍MH算法(Metropolis Hastings)
相关参考:
随机矩阵 - 百度百科
机器学习-白板推导系列-蒙特卡洛方法6 (Monte Carlo Method)- 平稳分布
机器学习-白板推导系列-蒙特卡洛方法6 (Monte Carlo Method)- 马尔可夫链(Markov Chain)