Loading [MathJax]/jax/element/mml/optable/SuppMathOperators.js

无记忆的随机过程——马尔科夫链的Python分析

马尔科夫链是一种用于描述系统从一个状态转移到另一个状态的随机过程。马尔科夫链的一个关键特性是无记忆性,即未来状态的概率只依赖于当前状态,而不依赖于过去的状态,这种性质使得马尔科夫链在许多领域中具有广泛的应用。马尔科夫链的概念由安德雷·马尔科夫在1906年提出,最初用于处理文学作品中的字符序列的统计特性。马尔科夫链的理论基础在20世纪得到了显著的发展,随着计算技术的进步,马尔科夫链的计算方法得到了极大的提升,进一步推动了其应用范围的扩展。20世纪下半叶,马尔科夫链在经济学、金融学和统计学中的应用逐渐成熟。尤其是随着贝叶斯统计和机器学习的兴起,隐马尔科夫模型(HMM)和马尔科夫链蒙特卡罗方法(MCMC)成为了数据分析和建模的重要工具。在现代,马尔科夫链的应用已渗透到各个科学和工程领域,并继续通过新的算法和计算能力的提升,不断扩展其应用范围。

一、马尔科夫链概述

在机器学习算法中,马尔可夫链(Markov chain)又称为离散时间马尔可夫链(discrete-time Markov chain),为状态空间中经过从一个状态到另一个状态的转换的随机过程。该过程要求具备“无记忆”的性质:下一状态的概率分布只能由当前状态决定,在时间序列中它前面的事件均与之无关。在马尔可夫链的每一步,系统根据概率分布,可以从一个状态变到另一个状态,也可以保持当前状态。状态的改变叫做转移,与不同的状态改变相关的概率叫做转移概率。随机漫步就是马尔可夫链的例子,随机漫步中每一步的状态是在图形中的点,每一步可以移动到任何一个相邻的点。

1.1 马尔科夫链的引入

搜索引擎采用不同的算法将网页按给定的关键字以相关度递减的方式排序.其中一种算法的思想是计算马尔可夫链的稳态分布。互联网由一系列相互链接的网页组成。这些网页和它们之间的链接关系构成了一张图。如图1所示,图中的节点是所有的网页X 。若网页i有一个到网页j的链接,则图中有一条由ij的弧(有向边)。

图1 图2

在图1中 P(A,B)=1/2P(D,E)=1/3可由这些网页节点A、B、C、D和E的出度来导出,如A节点有两个出度为2,一条去B,一条去D,所以P(A,B)=1/2。直观上看,一个高级别网页所指向的网页也应具有较高的级别。因此,我们定义网页 i级别 π(i) 为一个正数,并且满足以下公式:

π(i)=jXπ(j)P(j,i),iX

P(j,i) 表示所有网页j的外向链接中指向i的链接所占的比例。如果j没有指向i的链接,则P(j,i)=0。可以将这些等式以矩阵的形式记作:

π=πP

这里,π 是一个以π(i)为分量的行向量,而 P 是一个以P(i,j)为元素的方阵。式(1)称为平衡方程。这里我们注意到,如果一个向量π是这个方程的解,那么π 的任意倍数也是这个方程的解。为方便起见,我们将解归一化,使得所有网页的级别之和为1:

iXπ(i)=1

对于图1的例子,平衡方程为:

π(A)=π(C)+π(D)×(1/3)π(B)=π(A)×(1/2)+π(D)×(1/3)+π(E)×(1/2)π(C)=π(B)+π(E)×(1/2)π(D)=π(A)×(1/2)π(E)=π(D)×(1/3).

再加上π(i)的和为1这一条件,可以得到:

π=[π(A),π(B),π(C),π(D),π(E)]=139[12,9,10,6,2]

由此可以看到,网页A 具有最高的级别,网页E具有最低的级别。运用这一方法的搜索引擎会将这些网页级别与其他因素相结合来进行网页排序。一些搜索引擎也会采用这一度量的变种来对网页进行排序,这就是谷歌公司创始人之一拉里·佩奇提出的PageRank排序原理。

【平衡方程描述了系统在稳态分布下的流入和流出状态的平衡关系。具体而言,对于每个状态i,平衡方程要求系统进入状态i 的总概率等于系统从状态 i 转移出去的总概率。设πi 是状态 i 的稳态概率,Pji是从状态j 转移到状态 i 的概率,则平衡方程为:

πi=jπjPji

这意味着,对于每个状态i,从所有其他状态j转移到状态 i 的概率总和等于在稳态分布下状态 i 的概率πi​。同时,这也意味着在稳态分布下,系统在任意状态的概率保持不变。
假设有一个三状态的马尔科夫链,其转移矩阵P 为:

P=(0.20.30.50.10.60.30.40.50.1)

我们希望找到稳态分布 π=(π1,π2,π3),使得 πP=π。平衡方程可以表示为:

{π1=0.2π1+0.1π2+0.4π3π2=0.3π1+0.6π2+0.5π3π3=0.5π1+0.3π2+0.1π3

同时,我们需要满足概率归一化条件:

π1+π2+π3=1

通过求解这些方程,可以得到系统的稳态分布。这种分布描述了系统在长期运行后各个状态的概率,并且这些概率不会随着时间的推移而改变。】

1.2 马尔科夫链的理论

概率向量: 若向量 u=(u1u2un)T 满足 ui\sum_{i=1}^n u_i=1, 则称 \boldsymbol{u} 为概率向量。若矩阵 \left[a_{i j}\right]_{n \times n} 各个行向量均为概率向量, 称矩阵 \left[a_{i j}\right]_{n \times n} 为随机(概率)矩阵。
三阶方阵 \boldsymbol{P}=\left[\begin{array}{ccc}0.8 & 0.15 & 0.05 \\ 0.2 & 0.5 & 0.3 \\ 0.1 & 0.2 & 0.7\end{array}\right] 就是概率矩阵。
马尔科夫链:X_1, X_2, \cdots, X_n, \cdots 为离散的随机变量序列, 简记为 \left\{X_n\right\}, X_n 的所有可能取值的全体称为 \left\{X_n\right\} 的状态空间, 记为 E=\left\{x_1, x_2, x_3, \cdots\right\} 。若对任意的正整数 n 及任意的 x_{i_1}, x_{i_1}, x_{i_1}, \cdots, x_{i_n}, x_{i_{n+1}} \in E, 只要

P\left(X_1=x_{i_1}, X_2=x_{i_2}, X_3=x_{i_3}, \cdots, X_n=x_{i_n}\right)>0

就有

P\left(X_{n+1}=x_{i_{e+1}} \mid X_1=x_{i_i}, X_2=x_{i_2}, X_3=x_{i_1}, \cdots, X_n=x_{i_n}\right)=P\left(X_{n+1}=x_{i_{n+1}} \mid X_n=x_{i_e}\right)

则称 \left\{X_n\right\} 为马尔科夫链, 简称马氏链。
齐次马氏链: \left\{X_n\right\} 为马氏链, 若对任意的 x_i, x_j \in E, 总有

P\left(X_{n+1}=x_j \mid X_n=x_i\right)=P\left(X_{m+1}=x_j \mid X_m=x_i\right)

则称 \left\{X_n\right\} 为齐次马氏链。
转移矩阵:\left\{X_n\right\} 为齐次马氏链, 称 P\left(X_{n+k}=x_j \mid X_n=x_i\right)\left\{X_n\right\} 从状态 x_i 到状态 x_jk 步转移概率, 记作 p_{i j}(k); 称以 p_{i j}(k)\left(x_i, x_j \in E\right) 为元素的矩阵为 \left\{X_n\right\}k 步转移矩阵, 记作 \boldsymbol{P}(k), 特别地, 将一步转移概率和一步转移矩阵分别记为 p_{i j}\boldsymbol{P} 。显然, 一步和多步转移矩阵都是随机矩阵。
极限矩阵:若矩阵 \boldsymbol{A}(n)=\left[\begin{array}{cccc}a_{11}(n) & a_{12}(n) & \cdots & a_{1 m}(n) \\ a_{21}(n) & a_{22}(n) & \cdots & a_{2 m}(n) \\ \cdots & \cdots & \cdots & \cdots \\ a_{m 1}(n) & a_{m 2}(n) & \cdots & a_{m m}(n)\end{array}\right] 每一个元素 a_{i j}(n) (n=1,2,3, \cdots) 都是数列 \left\{a_{i j}(n)\right\} 的项, 称矩阵 \boldsymbol{A}(n) 为数列矩阵; 若对任意的 i, j=1,2, \cdots, m 都有 \lim _{n \rightarrow \infty} a_{i j}(n)=a_{i j} 成立, 即每个数列的极限都存在, 称矩阵 \boldsymbol{A}(n)n \rightarrow \infty 时的极限为 \boldsymbol{A}=\left[a_{i j}\right], \boldsymbol{A}\boldsymbol{A}(n) 的极限矩阵。
平稳分布(或稳态分布):是指在长期运行下,马尔科夫链中各状态的概率分布保持不变。若一个状态从某个分布开始,经过多次转移后,最终达到一个稳定的概率分布,此分布即为平稳分布。数学上,平稳分布\pi满足前面提到的平衡方程(1):\pi P = \pi,其中,\pi 是行向量,P 是状态转移矩阵。

齐次马氏链的 k 步转移矩阵 \boldsymbol{P}(k)=\boldsymbol{P}^k
若齐次马氏链的 k 步转移矩阵 \boldsymbol{P}(k)k \rightarrow \infty 的极限矩阵存在, 记为 \lim \boldsymbol{P}, 则 \lim \boldsymbol{P} 是一个随机矩阵。该随机矩阵也是一个转移矩阵, 且对任意的自然数 n, 都有 (\lim \boldsymbol{P})^n=\lim \boldsymbol{P}
若令P^*=\lim \boldsymbol{P}为极限矩阵,即

P^* = \lim_{k \to \infty} P^k

P^* 满足:

P∗P^* = P^**P=P^*

这意味着从P^*的每一行都是相同的平稳分布\pi
若齐次马氏链的 k 步转移矩阵 \boldsymbol{P}(k) 的极限矩阵存在, 则随着系统的不断演化, 最终系统状态之间的转移概率保持不变, 系统体现出统计规律性的特征, 就演化为一个稳定的系统。为了研究问题方便起见, 在实际问题中经常考虑有限介状态的齐次马氏链, 它们的 k 步转移矩阵的极限矩阵总存在。
假设抽烟和不抽烟的人群之间有一定的概率互相转化,且相互转化的概率如下图a所示。抽烟和不抽烟初始人数[500,500],则抽烟人数随着天数的变化见下图b。

图a 图b
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import math
transfer_matrix = np.array([[0.8, 0.1],[0.2,0.9]],dtype='float32')
somking_matrix = np.array([[500],[500]],dtype='float32')
number_somking = [500]     #初始人数分布
number_nosomking = [500]

for i in range(30):
    somking_matrix = np.dot(transfer_matrix,somking_matrix)
    #print(somking_matrix)
    number_somking.append(somking_matrix[0][0])
    number_nosomking.append(somking_matrix[1][0])
    
print(number_somking)
print(number_nosomking)

x = np.arange(31)
plt.plot(x,number_somking,label='somking')
plt.plot(x,number_nosomking,label='nosomking')
plt.legend()
plt.show()

二、马尔科夫链实例

2.1 案例1

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

图3 图4
import numpy as np
import matplotlib.pyplot as plt

# 状态转移矩阵
transfer_matrix = np.array([[0.9, 0.075, 0.025],
                             [0.15, 0.8, 0.05],
                             [0.25, 0.25, 0.5]])

# 状态名和对应的索引位置
states = ["Bull", "Bear", "Stagnant"]
positions = np.array([1, 2, 3])  # 状态的图形位置

# 计算特征向量
eigenvalues, eigenvectors = np.linalg.eig(transfer_matrix)
max_eigval_index = np.argmax(np.abs(eigenvalues))
steady_state_vector = eigenvectors[:, max_eigval_index].real

# 归一化得到极限分布
steady_state_distribution = steady_state_vector / steady_state_vector.sum()

# 打印状态的极限分布,保留两位小数
print("状态的极限分布(保留两位小数):")
for state, value in zip(states, np.around(steady_state_distribution, 2)):
    print(f"{state}: {value:.2f}")

2.2 案例2

保险公式将人的生存状态定义为三种:健康、疾病和死亡(第三中状态),用Xn=3表示.今年健康,明年可能突发疾病或偶然事件而死亡;今年疾病,明年可能继续疾病、健康或死亡,而一旦死亡,则不能再转为健康或疾病。根据统计,三种状态相互转化的概率如图5所示。

图5 图6

2.3 案例3

图7 图8

三、Python实现

#上面转移矩阵的极限分布计算
import numpy as np

# 定义转移概率矩阵
P = np.array([[0.2, 0.3, 0.5],
              [0.1, 0.6, 0.3],
              [0.4, 0.5, 0.1]])

# 计算多步转移矩阵
steps = [1, 2, 5, 10, 50, 100]
for step in steps:
    P_n = np.linalg.matrix_power(P, step)
    print(f"P^{step}:\n{np.round(P_n, 2)}\n")

# 计算极限矩阵
eigvals, eigvecs = np.linalg.eig(P.T)
eigvec1 = eigvecs[:, np.isclose(eigvals, 1)]

# 将特征向量归一化以使其总和为1
stationary = eigvec1 / eigvec1.sum()
stationary = stationary.real.flatten()
stationary = np.round(stationary, 2)
print("Stationary distribution:\n", stationary)
P^2:
[[0.27 0.49 0.24]
 [0.2  0.54 0.26]
 [0.17 0.47 0.36]]

P^10:
[[0.21 0.51 0.28]
 [0.21 0.51 0.28]
 [0.21 0.51 0.28]]

P^50:
[[0.21 0.51 0.28]
 [0.21 0.51 0.28]
 [0.21 0.51 0.28]]

Stationary distribution:
 [0.21 0.51 0.28]

总结

马尔科夫链作为一种描述随机过程的数学工具,在多个领域中具有广泛的应用。其核心概念是无记忆性,使得它可以简洁而有效地建模复杂系统的动态行为。通过转移矩阵和极限矩阵,马尔科夫链能够描述系统在长期运行中的稳定状态。随着计算技术的进步,马尔科夫链的计算和模拟方法得到了极大的提升,进一步推动了其在经济学、生物信息学、通信、计算机科学以及物理化学等领域的应用和发展。马尔科夫链不仅在理论研究中具有重要地位,而且在实际问题解决中也展现出强大的实用性和灵活性。未来,随着数据规模的不断增长和计算能力的提升,马尔科夫链的应用将会更加广泛和深入,继续为科学研究和工程实践提供强有力的支持。

参考文献

  1. 第1章 PageRank—A
  2. 马尔可夫链预测 (Markov Chain)
  3. MarcovChain介绍
  4. 马尔可夫链预测 (Markov Chain)
posted @   郝hai  阅读(487)  评论(0编辑  收藏  举报
相关博文:
阅读排行:
· DeepSeek智能编程
· 精选4款基于.NET开源、功能强大的通讯调试工具
· [翻译] 为什么 Tracebit 用 C# 开发
· 腾讯ima接入deepseek-r1,借用别人脑子用用成真了~
· DeepSeek崛起:程序员“饭碗”被抢,还是职业进化新起点?
点击右上角即可分享
微信分享提示