最簡單的隨機過程——馬爾科夫鏈的Python分析

郝hai發表於2024-06-15

馬爾科夫鏈是一種用於描述系統從一個狀態轉移到另一個狀態的隨機過程。它得名於俄羅斯數學家安德雷·馬爾科夫,他在20世紀初提出了這種數學模型。馬爾科夫鏈的一個關鍵特性是無記憶性,即未來狀態的機率只依賴於當前狀態,而不依賴於過去的狀態。這種性質使得馬爾科夫鏈在許多領域中具有廣泛的應用,如統計學、物理學、生物學、經濟學、電腦科學等。馬爾科夫鏈的概念由安德雷·馬爾科夫在1906年提出,最初用於處理文學作品中的字元序列的統計特性。馬爾科夫希望透過研究無記憶性的隨機過程,挑戰傳統的獨立事件統計學理論。馬爾科夫鏈的理論基礎在20世紀得到了顯著的發展,特別是在第二次世界大戰期間和戰後時期。隨著計算技術的進步,馬爾科夫鏈的計算方法得到了極大的提升,進一步推動了其應用範圍的擴充套件。20世紀下半葉,馬爾科夫鏈在經濟學、金融學和統計學中的應用逐漸成熟。尤其是隨著貝葉斯統計和機器學習的興起,隱馬爾科夫模型(HMM)和馬爾科夫鏈蒙特卡羅方法(MCMC)成為了資料分析和建模的重要工具。在現代,馬爾科夫鏈的應用已滲透到各個科學和工程領域,並繼續透過新的演算法和計算能力的提升,不斷擴充套件其應用範圍。例如,強化學習中的馬爾科夫決策過程(MDP)是人工智慧研究中的一個重要方向,用於設計和分析智慧代理在動態環境中的決策策略。

一、馬爾科夫鏈概述

在機器學習演算法中,馬爾可夫鏈(Markov chain)是個很重要的概念。馬爾可夫鏈(Markov chain),又稱離散時間馬爾可夫鏈(discrete-time Markov chain),為狀態空間中經過從一個狀態到另一個狀態的轉換的隨機過程。該過程要求具備“無記憶”的性質:下一狀態的機率分佈只能由當前狀態決定,在時間序列中它前面的事件均與之無關。這種特定型別的“無記憶性”稱作馬爾可夫性質。馬爾科夫鏈是指數學中具有馬爾科夫性質的離散事件隨機過程。在馬爾可夫鏈的每一步,系統根據機率分佈,可以從一個狀態變到另一個狀態,也可以保持當前狀態。狀態的改變叫做轉移,與不同的狀態改變相關的機率叫做轉移機率。隨機漫步就是馬爾可夫鏈的例子。隨機漫步中每一步的狀態是在圖形中的點,每一步可以移動到任何一個相鄰的點,在這裡移動到每一個點的機率都是相同的(無論之前漫步路徑是如何的)。

1.1 馬爾科夫鏈的引入

搜尋引擎採用不同的演算法將網頁按給定的關鍵字以相關度遞減的方式排序.其中一種演算法的思想是計算馬爾可夫鏈的穩態分佈。網際網路由一系列相互連結的網頁組成。這些網頁和它們之間的連結關係構成了一張圖。如圖1所示,圖中的節點是所有的網頁\(\mathcal{X}\) 。若網頁\(i\)有一個到網頁\(j\)的連結,則圖中有一條由\(i\)\(j\)的弧(有向邊)。

圖1 圖2

在圖1中 \(P(A, B) = 1/2 ,P(D, E) = 1/3\)可由這些網頁節點A、B、C、D和E的出度來匯出。直觀上看,一個高階別網頁所指向的網頁也應具有較高的級別。因此,我們定義網頁 \(i\)級別 $ \pi(i)$ 為一個正數,並且滿足以下公式:

\[\pi(i) = \sum_{j \in \mathcal{X}} \pi(j) P(j, i), \quad i \in \mathcal{X} \]

\(P(j, i)\) 表示所有網頁\(j\)的外向連結中指向\(i\)的連結所佔的比例。如果\(j\)沒有指向\(i\)的連結,則$P(j, i) = 0 $。可以將這些等式以矩陣的形式記作:

\[\pi = \pi P \tag{1} \]

這裡,$\pi $ 是一個以$\pi(i) $為分量的行向量,而 \(P\) 是一個以$P(i, j) $為元素的方陣。

式(1)稱為平衡方程。這裡我們注意到,如果一個向量\(\pi\)是這個方程的解,那麼\(\pi\) 的任意倍數也是這個方程的解。為方便起見,我們將解歸一化,使得所有網頁的級別之和為1:

\[\sum_{i \in \mathcal{X}} \pi(i) = 1 \]

對於圖1的例子,平衡方程為:

\[\begin{aligned} &\pi(A) = \pi(C) + \pi(D) \times (1/3) \\ &\pi(B) = \pi(A) \times (1/2) + \pi(D) \times (1/3) + \pi(E) \times (1/2) \\ &\pi(C) = \pi(B) + \pi(E) \times (1/2) \\ &\pi(D) = \pi(A) \times (1/2) \\ &\pi(E) = \pi(D) \times (1/3). \end{aligned} \]

再加上\(\pi(i)\)的和為1這一條件,可以得到:

\[\pi = [\pi(A), \pi(B), \pi(C), \pi(D), \pi(E)] = \frac{1}{39} [12, 9, 10, 6, 2] \]

由此可以看到,網頁\(A\) 具有最高的級別,網頁\(E\)具有最低的級別。運用這一方法的搜尋引擎會將這些網頁級別與其他因素相結合來進行網頁排序。一些搜尋引擎也會採用這一度量的變種來對網頁進行排序,這就是谷歌公司創始人之一拉里·佩奇提出的PageRank排序原理。

1.2 馬爾科夫鏈的理論

機率向量: 若向量 \(\boldsymbol{u}=\left(u_1 u_2 \cdots u_n\right)^{\mathrm{T}}\) 滿足 \(u_i \geqslant 0(i=1,2, \cdots, n)\)\(\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_j\)\(k\) 步轉移機率, 記作 \(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)\) 的極限矩陣。

齊次馬氏鏈的 \(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}\)
若齊次馬氏鏈的 \(k\) 步轉移矩陣 \(\boldsymbol{P}(k)\) 的極限矩陣存在, 則隨著系統的不斷演化, 最終系統狀態之間的轉移機率保持不變, 系統體現出統計規律性的特徵, 就演化為一個穩定的系統。為了研究問題方便起見, 在實際問題中經常考慮有限介狀態的齊次馬氏鏈, 它們的 \(k\) 步轉移矩陣的極限矩陣總存在。

二、馬爾科夫鏈例項

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]

總結

馬爾科夫鏈作為一種描述隨機過程的數學工具,在多個領域中具有廣泛的應用。其核心概念是無記憶性,使得它可以簡潔而有效地建模複雜系統的動態行為。透過轉移矩陣和極限矩陣,馬爾科夫鏈能夠描述系統在長期執行中的穩定狀態。隨著計算技術的進步,馬爾科夫鏈的計算和模擬方法得到了極大的提升,進一步推動了其在經濟學、生物資訊學、通訊、電腦科學以及物理化學等領域的應用和發展。馬爾科夫鏈不僅在理論研究中具有重要地位,而且在實際問題解決中也展現出強大的實用性和靈活性。未來,隨著資料規模的不斷增長和計算能力的提升,馬爾科夫鏈的應用將會更加廣泛和深入,繼續為科學研究和工程實踐提供強有力的支援。

最簡單的隨機過程——馬爾科夫鏈的Python分析

參考文獻

  1. 第1章 PageRank—A
  2. 馬爾可夫鏈預測 (Markov Chain)
  3. 馬爾科夫決策
  4. 馬爾可夫鏈預測 (Markov Chain)

相關文章