香港中文大學數學系
2026年6月27日
# 四室老鼠馬爾可夫鏈 (Four-Compartment Mouse Markov Chain)
# 描述：
# 1. 定義轉移矩陣
# 2. 驗證隨機矩陣的性質
# 3. 計算隨時間變化的分佈
# 4. 尋找 3 步後在第 2 室的概率
# 5. 計算平穩分佈
# 6. 驗證平穩方程
import numpy as np
from scipy.linalg import eig
# ============================================================
# 轉移矩陣 (Transition Matrix)
# ============================================================
P = np.array([
[0, 2/3, 0, 1/3],
[2/3, 0, 1/3, 0 ],
[0, 1/2, 0, 1/2],
[1/2, 0, 1/2, 0 ]
], dtype=float)
# ============================================================
# 驗證矩陣 (Validate Matrix)
# ============================================================
def validate_transition_matrix(P):
"""
檢查矩陣是否為有效的轉移矩陣。
"""
rowsum = P.sum(axis=1)
if not np.allclose(rowsum, 1):
raise ValueError(
f"Invalid transition matrix. Row sums = {rowsum}"
)
if np.any(P < 0):
raise ValueError(
"Transition matrix contains negative probabilities."
)
print("✓ Transition matrix is valid.\n")
validate_transition_matrix(P)## ✓ Transition matrix is valid.# ============================================================
# 初始狀態 (Initial State)
# ============================================================
# 從第 2 室開始
pi0 = np.array([0, 1, 0, 0])
print("Initial distribution:")## Initial distribution:## [0 1 0 0]# ============================================================
# 計算分佈 (Compute Distributions)
# ============================================================
def propagate_distribution(initial_dist, P, n_steps):
"""
計算 n_steps 步後的分佈。
π(n) = π(0) P^n
"""
distributions = [initial_dist]
current = initial_dist.copy()
for _ in range(n_steps):
current = current @ P
distributions.append(current)
return distributions
distributions = propagate_distribution(pi0, P, 3)
for step, dist in enumerate(distributions):
print(f"Distribution at step {step}:")
print(dist)
print()## Distribution at step 0:
## [0 1 0 0]
##
## Distribution at step 1:
## [0.66666667 0. 0.33333333 0. ]
##
## Distribution at step 2:
## [0. 0.61111111 0. 0.38888889]
##
## Distribution at step 3:
## [0.60185185 0. 0.39814815 0. ]# ============================================================
# 3 步後狀態 2 的概率 (Probability of state 2 after 3 steps)
# ============================================================
pi3 = distributions[3]
prob_state_2 = pi3[1]
print("=" * 60)## ============================================================## Probability mouse is in compartment 2 after 3 steps## ============================================================## 0.0# ============================================================
# 直接計算 P^3 (Compute P^3 directly)
# ============================================================
P3 = np.linalg.matrix_power(P, 3)
print("=" * 60)## ============================================================## P^3## ============================================================## [[0. 0.60185185 0. 0.39814815]
## [0.60185185 0. 0.39814815 0. ]
## [0. 0.59722222 0. 0.40277778]
## [0.59722222 0. 0.40277778 0. ]]# ============================================================
# 平穩分佈 (Stationary Distribution)
# ============================================================
def stationary_distribution(P):
"""
解出 πP = π
使用與特徵值 1 相關聯的特徵向量。
"""
eigenvalues, eigenvectors = eig(P.T)
idx = np.argmin(np.abs(eigenvalues - 1))
stationary = np.real(eigenvectors[:, idx])
stationary = stationary / stationary.sum()
stationary = stationary.real
return stationary
pi_star = stationary_distribution(P)
print("=" * 60)## ============================================================## Stationary Distribution## ============================================================## [0.3 0.3 0.2 0.2]# ============================================================
# 驗證 (Verification)
# ============================================================
verification = pi_star @ P
print("π*P:")## π*P:## [0.3 0.3 0.2 0.2]## Difference:## [ 1.66533454e-16 -5.55111512e-17 -5.55111512e-17 -8.32667268e-17]if np.allclose(verification, pi_star):
print("✓ Stationary distribution verified.")
else:
print("✗ Verification failed.")## ✓ Stationary distribution verified.# ============================================================
# 第 2 室的長期概率 (Long-run probability of compartment 2)
# ============================================================
print("=" * 60)## ============================================================## Long-run probability of compartment 2## ============================================================## 0.3# 計算：
# P(X1=2, X4=1, X6=1, X18=1 | X0=1)
# 使用矩陣冪。
# 轉移矩陣：
# 目的地
# 1 2
# ----------------
# 1 | 0.5 0.5
# 2 | 0.3 0.7
import numpy as np
# =============================================================
# 轉移矩陣 (Transition Matrix)
# =============================================================
P = np.array(
[
[0.5, 0.5],
[0.3, 0.7]
],
dtype=float
)
# =============================================================
# 驗證函數 (Validation Functions)
# =============================================================
def validate_transition_matrix(matrix):
"""
驗證馬爾可夫轉移矩陣。
條件：
1. 必須是方陣
2. 無負概率
3. 每一列的總和為 1
"""
rows, cols = matrix.shape
if rows != cols:
raise ValueError(
"Transition matrix must be square."
)
if np.any(matrix < 0):
raise ValueError(
"Transition matrix contains negative entries."
)
row_sums = matrix.sum(axis=1)
if not np.allclose(row_sums, 1.0):
raise ValueError(
f"Rows do not sum to 1.\nRow sums = {row_sums}"
)
print("✓ Transition matrix is valid.\n")
# =============================================================
# 計算平穩分佈 (Compute Stationary Distribution)
# =============================================================
def stationary_distribution(matrix):
"""
使用特徵向量求解：
πP = π
"""
eigenvalues, eigenvectors = np.linalg.eig(
matrix.T
)
idx = np.argmin(
np.abs(eigenvalues - 1)
)
stationary = np.real(
eigenvectors[:, idx]
)
stationary = stationary / stationary.sum()
return stationary
# =============================================================
# 主程式 (Main Program)
# =============================================================
def main():
print("=" * 60)
print("TWO-STATE MARKOV CHAIN")
print("=" * 60)
print()
# ---------------------------------------------------------
# 驗證矩陣
# ---------------------------------------------------------
validate_transition_matrix(P)
print("Transition Matrix P")
print(P)
print()
# ---------------------------------------------------------
# 矩陣冪
# ---------------------------------------------------------
P2 = np.linalg.matrix_power(P, 2)
P3 = np.linalg.matrix_power(P, 3)
P12 = np.linalg.matrix_power(P, 12)
print("=" * 60)
print("P^2")
print("=" * 60)
print(P2)
print()
print("=" * 60)
print("P^3")
print("=" * 60)
print(P3)
print()
print("=" * 60)
print("P^12")
print("=" * 60)
print(P12)
print()
# ---------------------------------------------------------
# 提取所需的概率
# ---------------------------------------------------------
# P(X1=2 | X0=1)
p12 = P[0, 1]
# P(X4=1 | X1=2)
p3_21 = P3[1, 0]
# P(X6=1 | X4=1)
p2_11 = P2[0, 0]
# P(X18=1 | X6=1)
p12_11 = P12[0, 0]
print("=" * 60)
print("Required Entries")
print("=" * 60)
print(f"P12 = {p12:.12f}")
print(f"(P^3)21 = {p3_21:.12f}")
print(f"(P^2)11 = {p2_11:.12f}")
print(f"(P^12)11 = {p12_11:.12f}")
print()
# ---------------------------------------------------------
# 最終聯合概率
# ---------------------------------------------------------
joint_probability = (
p12
* p3_21
* p2_11
* p12_11
)
print("=" * 60)
print("FINAL JOINT PROBABILITY")
print("=" * 60)
print(
"P(X1=2, X4=1, X6=1, X18=1 | X0=1)"
)
print(f"= {joint_probability:.15f}")
print()
# ---------------------------------------------------------
# 平穩分佈
# ---------------------------------------------------------
pi_star = stationary_distribution(P)
print("=" * 60)
print("STATIONARY DISTRIBUTION")
print("=" * 60)
print(pi_star)
print()
print("Verification πP:")
print(pi_star @ P)
print()
if np.allclose(pi_star @ P, pi_star):
print("✓ Verified: πP = π")
else:
print("✗ Verification failed")
print()
# =============================================================
# 執行程式
# =============================================================
if __name__ == "__main__":
main()## ============================================================
## TWO-STATE MARKOV CHAIN
## ============================================================
##
## ✓ Transition matrix is valid.
##
## Transition Matrix P
## [[0.5 0.5]
## [0.3 0.7]]
##
## ============================================================
## P^2
## ============================================================
## [[0.4 0.6 ]
## [0.36 0.64]]
##
## ============================================================
## P^3
## ============================================================
## [[0.38 0.62 ]
## [0.372 0.628]]
##
## ============================================================
## P^12
## ============================================================
## [[0.375 0.625]
## [0.375 0.625]]
##
## ============================================================
## Required Entries
## ============================================================
## P12 = 0.500000000000
## (P^3)21 = 0.372000000000
## (P^2)11 = 0.400000000000
## (P^12)11 = 0.375000002560
##
## ============================================================
## FINAL JOINT PROBABILITY
## ============================================================
## P(X1=2, X4=1, X6=1, X18=1 | X0=1)
## = 0.027900000190464
##
## ============================================================
## STATIONARY DISTRIBUTION
## ============================================================
## [0.375 0.625]
##
## Verification πP:
## [0.375 0.625]
##
## ✓ Verified: πP = πimport numpy as np
# -------------------------------
# 步驟 1：定義轉移矩陣
# 行 (Rows) = 當前狀態，列 (Columns) = 下一個狀態
# 狀態：[多雲 (Cloudy), 晴天 (Sunny), 下雪 (Snowy)]
# -------------------------------
transmat = np.array([
[0.57, 0.22, 0.21], # 多雲 (Cloudy) 列
[0.21, 0.44, 0.10], # 晴天 (Sunny) 列
[0.21, 0.33, 0.33] # 下雪 (Snowy) 列
])
# -------------------------------
# 步驟 2：定義初始狀態向量
# 2020年11月12日 = 多雲
# -------------------------------
state0 = np.array([1, 0, 0]) # 多雲 = 1, 晴天 = 0, 下雪 = 0
# -------------------------------
# 步驟 3：預測函數
# state_n = transmat^n * state0
# -------------------------------
def forecast(transmat, state0, n_steps):
# 將轉移矩陣提升到 n 次方
transmat_n = np.linalg.matrix_power(transmat, n_steps)
# 乘以初始狀態向量
state_n = transmat_n @ state0
return transmat_n, state_n
# -------------------------------
# 步驟 4：計算未來 1、2、3 天的預測
# -------------------------------
for n in range(1, 4):
transmat_n, state_n = forecast(transmat, state0, n)
print(f"\nTransition Matrix ^ {n}:\n", transmat_n)
print(f"State vector at day {n}:\n", state_n)##
## Transition Matrix ^ 1:
## [[0.57 0.22 0.21]
## [0.21 0.44 0.1 ]
## [0.21 0.33 0.33]]
## State vector at day 1:
## [0.57 0.21 0.21]
##
## Transition Matrix ^ 2:
## [[0.4152 0.2915 0.211 ]
## [0.2331 0.2728 0.1211]
## [0.2583 0.3003 0.186 ]]
## State vector at day 2:
## [0.4152 0.2331 0.2583]
##
## Transition Matrix ^ 3:
## [[0.342189 0.289234 0.185972]
## [0.215586 0.211277 0.116194]
## [0.249354 0.250338 0.145653]]
## State vector at day 3:
## [0.342189 0.215586 0.249354]# -------------------------------
# 步驟 5：結果解釋
# -------------------------------
labels = ["Cloudy", "Sunny", "Snowy"]
transmat3, state3 = forecast(transmat, state0, 3)
print("\n--- Forecast for November 15 (3 days ahead) ---")##
## --- Forecast for November 15 (3 days ahead) ---## Cloudy: 0.34 (34%)
## Sunny: 0.22 (22%)
## Snowy: 0.25 (25%)# 來自表格的轉移矩陣（牛市 (Bull)、熊市 (Bear)、停滯 (Stagnant)）
import numpy as np
# -------------------------------
# 步驟 1：定義轉移矩陣
# 行 (Rows) = 當前狀態，列 (Columns) = 下一個狀態
# 狀態：[牛市, 熊市, 停滯]
# -------------------------------
transmat = np.array([
[0.900, 0.075, 0.025], # 牛市 (Bull) 列
[0.150, 0.800, 0.050], # 熊市 (Bear) 列
[0.250, 0.250, 0.500] # 停滯 (Stagnant) 列
])
# -------------------------------
# 步驟 2：定義初始狀態行向量
# 本週 = 熊市
# -------------------------------
state0 = np.array([0, 1, 0]) # 牛市=0, 熊市=1, 停滯=0
# -------------------------------
# 步驟 3：預測函數
# state_n = state0 * transmat^n
# -------------------------------
def forecast(transmat, state0, n_steps):
# 將轉移矩陣提升到 n 次方
transmat_n = np.linalg.matrix_power(transmat, n_steps)
# 將行向量乘以矩陣冪
state_n = state0 @ transmat_n
return transmat_n, state_n
# -------------------------------
# 步驟 4：計算未來 1、5、52、99 週的預測
# -------------------------------
for n in [1, 5, 52, 99]:
transmat_n, state_n = forecast(transmat, state0, n)
print(f"\nTransition Matrix ^ {n}:\n", transmat_n)
print(f"State vector at week {n}:\n", state_n)##
## Transition Matrix ^ 1:
## [[0.9 0.075 0.025]
## [0.15 0.8 0.05 ]
## [0.25 0.25 0.5 ]]
## State vector at week 1:
## [0.15 0.8 0.05]
##
## Transition Matrix ^ 5:
## [[0.70683 0.238305 0.054865]
## [0.47661 0.450515 0.072875]
## [0.54865 0.364375 0.086975]]
## State vector at week 5:
## [0.47661 0.450515 0.072875]
##
## Transition Matrix ^ 52:
## [[0.62500006 0.31249994 0.06249999]
## [0.62499988 0.31250011 0.06250001]
## [0.62499995 0.31250005 0.0625 ]]
## State vector at week 52:
## [0.62499988 0.31250011 0.06250001]
##
## Transition Matrix ^ 99:
## [[0.625 0.3125 0.0625]
## [0.625 0.3125 0.0625]
## [0.625 0.3125 0.0625]]
## State vector at week 99:
## [0.625 0.3125 0.0625]# -------------------------------
# 步驟 5：計算平穩分佈
# 解出 π 使得 π * transmat = π 且 sum(π)=1
# -------------------------------
# 轉置以求解線性方程組：(transmat.T - I)π = 0
A = transmat.T - np.eye(3)
# 加入歸一化約束條件
A = np.vstack([A, np.ones(3)])
b = np.array([0, 0, 0, 1])
# 求解最小二乘法系統
pi, _, _, _ = np.linalg.lstsq(A, b, rcond=None)
print("\n--- Stationary Distribution ---")##
## --- Stationary Distribution ---labels = ["Bull", "Bear", "Stagnant"]
for label, prob in zip(labels, pi):
print(f"{label}: {prob:.2f} ({prob*100:.0f}%)")## Bull: 0.63 (63%)
## Bear: 0.31 (31%)
## Stagnant: 0.06 (6%)