一个简单的马尔可夫过程例子

最新推荐文章于 2025-10-14 20:07:52 发布

原创最新推荐文章于 2025-10-14 20:07:52 发布 · 3.9w 阅读

87 ·

CC 4.0 BY-SA版权

Matlab 专栏收录该内容

4 篇文章

订阅专栏

本文介绍了马尔可夫过程的基本概念，包括随机过程、马尔可夫过程、马尔可夫链及其转移矩阵，并通过实例展示了如何求解稳态分布。

部署运行你感兴趣的模型镜像

什么是马尔可夫过程(Markov Process)

要说什么是马尔可夫过程，首先必须讲讲什么是随机过程(Stochastic Process)。

設 $(\Omega, \mathcal{F}, P)$ 為一概率空間，另設集合 $T$ 為一指標集合。如果對於所有 $t\in T$ ，均有一隨機變量 $\xi_t(\omega)$ 定義於概率空間 $(\Omega, \mathcal{F}, P)$ ，則集合 $\{\xi_t(\omega)|t\in T\}$ 為一隨機過程。—Wikipedia

其实简单来说，随机过程就是对一连串随机变量(或事件)变迁或者说动态关系的描述。

马尔可夫过程其实就是随机过程的一种：

A stochastic process ${s_t:t\in N}$ is said to be Markov if

p r o b (s t + 1 | s t) = p r o b (s t + 1 | s t)

$prob(s_{t+1}|s^t)=prob(s_{t+1}|s_t)$
where

st $s^t$ is the history up to period

t $t$ and

st $s_t$ is the realization of the state at period

t $t$ .

换句话说，马尔可夫过程说的就是在随机过程中，明天的状态只取决于今天的状态而和今天之前的状态无关。也就是不管你是怎样达到 $s_t$ 的，只要你今天的状态是 $s_t$ , 你明天状态的可能性就决定好了。

我们经常还会看到马尔可夫链。其实马尔可夫链就是状态空间为可数集的马尔可夫过程：

Let $s_t\in S$ and if $S$ is finite (countable), we call the Markov Process with this countable state space Markov Chain.

转移矩阵(Transition Matrix)

既然当知道 $s_t$ 时我们就能知道 $s_{t+1}$ 发生的概率，那么我们自然可以把有限可数的状态空间的转移法(transistion law)写成一个矩阵。举个例子，如果状态空间只包含了两个状态A和B。如果今天是A,那么明天还是A的概率为0.8。如果今日是B,那么明天还是B的概率为0.7。转移矩阵可写为：

P = ∣ ∣ ∣ 0.8 0.3 0.2 0.7 ∣ ∣ ∣

$\mathbf{P} = \begin{vmatrix} 0.8 & 0.2 \\\ 0.3 & 0.7 \end{vmatrix}$

稳态(steady－state)

首先我们先定义一个横向量 $\pi$ ；

π ＝ ∣ ∣ p r o b (A) p r o b (B) ∣ ∣

$\pi＝\begin{vmatrix} prob(A) & prob(B) \end{vmatrix}$

假设我们今天的状态是A而不是B,那么我们明天是A还是B的可能性为：

∣ ∣ 10 ∣ ∣ \times ∣ ∣ ∣ 0.8 0.3 0.2 0.7 ∣ ∣ ∣ = ∣ ∣ 0.8 0.2 ∣ ∣

$\begin{vmatrix} 1 & 0 \end{vmatrix} \times \begin{vmatrix} 0.8 & 0.2 \\\ 0.3 & 0.7 \end{vmatrix} = \begin{vmatrix} 0.8 & 0.2 \end{vmatrix}$

后天是A还是B的可能性为：

∣ ∣ 0.8 0.2 ∣ ∣ \times ∣ ∣ ∣ 0.8 0.3 0.2 0.7 ∣ ∣ ∣ = ∣ ∣ 0.7 0.3 ∣ ∣

$\begin{vmatrix} 0.8 & 0.2 \end{vmatrix} \times \begin{vmatrix} 0.8 & 0.2 \\\ 0.3 & 0.7 \end{vmatrix} = \begin{vmatrix} 0.7 & 0.3 \end{vmatrix}$

大后天是A还是B的可能性为：

∣ ∣ 0.7 0.3 ∣ ∣ \times ∣ ∣ ∣ 0.8 0.3 0.2 0.7 ∣ ∣ ∣ = ∣ ∣ 0.65 0.35 ∣ ∣

$\begin{vmatrix} 0.7 & 0.3 \end{vmatrix} \times \begin{vmatrix} 0.8 & 0.2 \\\ 0.3 & 0.7 \end{vmatrix} = \begin{vmatrix} 0.65 & 0.35 \end{vmatrix}$

那么会不会某一时刻后我们的横向量 $\pi$ 能到达一个不动点：

π * \times P = π *

$\pi^* \times \mathbf{P}=\pi^ *$
这么繁杂的工作还是要我们的小奴Matlab来完成吧。

%initial pi%
pi=[1,0];

%transistion matrix%
P=[0.8,0.2;0.3,0.7];

%critical value%
cv=1;

while cv>0.0000000001;
    pi2=pi*P;
    cv=abs(min(pi2-pi));
    pi=pi2;
end

得出的结果为:

π * = ∣ ∣ 0.6 0.4 ∣ ∣

$\pi^*=\begin{vmatrix} 0.6 & 0.4 \end{vmatrix}$

让我们换一个initial pi来试试看:

%initial pi%
pi=[0.35,0.65];

%transistion matrix%
P=[0.8,0.2;0.3,0.7];

%critical value%
cv=1;

while cv>10e-06;
    pi2=pi*P;
    cv=abs(min(pi2-pi));
    pi=pi2;
end

得出的结果仍然为:

π * = ∣ ∣ 0.6 0.4 ∣ ∣

$\pi^*=\begin{vmatrix} 0.6 & 0.4 \end{vmatrix}$

其实仔细观察一下就可以看出这个 $\pi$ 其实是转移矩阵 $\mathbf{P}$ 的转置的对应特征值(Eigenvalue)等于1的特征向量(Eigenvector)：

π = π P π (I - P) = 0 (I - P') π' = 0 (P' - I) π' = 0

$\begin{aligned} \\ & \pi = \pi \mathbf{P} \\\ & \pi(\mathbf{I}-\mathbf{P}) = 0 \\\ & (\mathbf{I}-\mathbf{P'})\pi' = 0 \\\ & (\mathbf{P'}-\mathbf{I})\pi' = 0 \end{aligned}$

既然知道 $\pi$ 其实是 $\mathbf{P‘}$ 的特征向量，我们就可以直接通过Matlab找特征向量而不必要用繁琐的迭代法啦。

%%%%%Eigenvector%%%%%
[n,v]=eig(P');
ppi=n(:,1);
scale=sum(ppi); %ensure the summation of the distribution vector pi is 1
ppi=ppi/scale;

得出的结果绝对是当然地为:

π * = ∣ ∣ 0.6 0.4 ∣ ∣

$\pi^*=\begin{vmatrix} 0.6 & 0.4 \end{vmatrix}$

一个简单的模型

假设我们在研究一个粒子的运动。这个粒子随机地在A门与B门之间跳动。如果这个粒子跳入A门，那么它下一次跳入A门的概率为0.8。如果这个粒子跳入B门，那么它下次跳入B门的概率为0.7。我们的问题是，现在有100000个这样的粒子让它们跳跃，让它们跳跃1000次,我们会观察到多少个粒子在A门多少个在B门？

思路：

首先随机生成100000个粒子。1代表该粒子在A门，0代表在B门。

 size=100000;
 a=round(rand(1,size));

生成10000个以概率0.8为1，概率0.2为0的随机数。同时在生成10000个以概率0.7为0，概率0.3为1随机数。

 AA=rand(1,size)<=0.8;
 BB=rand(1,size)<=0.3;

用2a－AA。如果值为1说明该粒子从A门跳转到A门；如果值为2则说明该粒子从A门跳转到B门。另外值为0，－1则说明该粒子初始状态是在B门。

用2a－BB。如果值为0说明该粒子从B门跳转到B门；如果值为-1则说明该粒子从B门跳转到A门。另外值为2，1则说明该粒子初始状态是在A门。

A=2*a-AA;
B=2*a-BB;

找到A中大于等于1的元素与B中小于等于0的元素。同时新建一个向量a2用于存储结果。

indxA=find(A>=1);
indxB=find(B<=0);

a2=zeros(1,size);

将A中大于等于1的元素与B中小于等于0的元素按照对应的位置放入新的向量a中。

for i=1:length(indxA);
    a2(1,indxA(1,i))=A(1,indxA(1,i));
end

for i=1:length(indxB);
    a2(1,indxB(1,i))=B(1,indxB(1,i));
end

这时a2中1和－1代表粒子在A门；2和0代表粒子在B门。将－1换成1，2换成0。

indxA=find(a2==-1);
indxB=find(a2==2);

for i=1:length(indxA);
    a2(1,indxA(1,i))=1;
end

for i=1:length(indxB);
    a2(1,indxB(1,i))=0;
end

通过while loop让粒子跳跃1000次。

%%%%%Example%%%%%

%generate observations% 
size=100000;
a=round(rand(1,size)); %each column is one particle%

j=1000;
while j>1
AA=rand(1,size)<=0.8;
BB=rand(1,size)<=0.3;

A=2*a-AA;
B=2*a-BB;

indxA=find(A>=1);
indxB=find(B<=0);

a2=zeros(1,size);

for i=1:length(indxA);
    a2(1,indxA(1,i))=A(1,indxA(1,i));
end

for i=1:length(indxB);
    a2(1,indxB(1,i))=B(1,indxB(1,i));
end

indxA=find(a2==-1);
indxB=find(a2==2);

for i=1:length(indxA);
    a2(1,indxA(1,i))=1;
end

for i=1:length(indxB);
    a2(1,indxB(1,i))=0;
end

a=a2;
j=j-1;
end

sum(a2)

得出的结果大概在60000左右。也就是说大概能观察到60000个粒子在A门，40000个在B门。

正如前面得出的 $\pi^*=\begin{vmatrix} 0.6 & 0.4 \end{vmatrix}$ ：粒子出现在A门的概率为0.6，出现在B门的概率为0.4。