马尔可夫链

概述

马尔可夫链（Markov Chain）是一种特殊的随机过程，具有无记忆性（Memorylessness）的特点——当前状态只与上一状态有关，与更早的历史无关。¹

马尔可夫链在搜索引擎（PageRank）、语音识别、自然语言处理、推荐系统、金融建模、统计物理等领域都有广泛应用。

---

马尔可夫性质

条件独立性

设 $\{X_t{\}}_{t\in \mathbb{T}}$ 为一个离散时间随机过程，若对任意时刻 $t$ 和任意状态序列 $i_0,i_1,\dots ,i_{t-1},i_t,j$，满足：

$$P(X_{t+1}=j\mid X_t=i_t,X_{t-1}=i_{t-1},\dots ,X_0=i_0)=P(X_{t+1}=j\mid X_t=i_t)$$

则称该过程具有马尔可夫性质（Markov Property）。

严格表述

更一般地，马尔可夫性质可以表述为：

$$P(X_{t+1}\in A\mid {\mathcal{F}}_t)=P(X_{t+1}\in A\mid X_t)$$

其中 ${\mathcal{F}}_t=\sigma (X_0,X_1,\dots ,X_t)$ 是到时刻 $t$ 为止的自然 $\sigma$-代数。

直观理解

未来只与现在有关，与过去无关。

这意味着，如果我们知道当前状态 $X_t$，那么关于过去的任何信息都不会帮助我们更好地预测 $X_{t+1}$。

---

离散时间马尔可夫链的严格定义

概率论框架

一个离散时间马尔可夫链（Discrete-Time Markov Chain, DTMC）是一个三元组 $(S,\mathcal{A},\mathbf{P})$，其中：

要素	符号	说明
状态空间	$S$	可数集合（有限或可数无限），包含所有可能状态
$\sigma$-代数	$\mathcal{A}$	$S$ 上的 Borel $\sigma$-代数 $\mathcal{B}(S)$
转移核	$\mathbf{P}$	满足 Kolmogorov 不等式的转移概率族

转移核的定义

转移核（Transition Kernel）是一个映射 $\mathbf{P}:S\times \mathcal{B}(S)\to [0,1]$，满足：

1. 对任意 $x\in S$，映射 $A\mapsto \mathbf{P}(x,A)$ 是 $(S,\mathcal{A})$ 上的概率测度
2. 对任意 $A\in \mathcal{A}$，映射 $x\mapsto \mathbf{P}(x,A)$ 是 $\mathcal{A}$-可测函数

对于离散状态空间 $S=\{s_1,s_2,\dots \}$，转移核退化为转移概率矩阵：

$$p_{ij}=\mathbf{P}(s_i,\{s_j\})=P(X_{t+1}=s_j\mid X_t=s_i)$$

转移矩阵的性质

转移概率矩阵 $\mathbf{P}$ 是一个 $n\times n$（若 $|S|=n$）或 $\infty \times \infty$（若 $|S|=\infty$）的行随机矩阵：

$$\mathbf{P}=\left(\begin{array}{cccc}{p}_{11}& p_{12}& \cdots & p_{1n}\\ p_{21}& p_{22}& \cdots & p_{2n}\\ ⋮& ⋮& \ddots & ⋮\\ p_{n1}& p_{n2}& \cdots & p_{nn}\end{array}\right)$$

其中 $p_{ij}\ge 0$ 且 ${\sum }_{j\in S}{p}_{ij}=1$。

初始分布

初始分布（Initial Distribution）$\mu =({\mu }_i{)}_{i\in S}$ 是一个概率向量：

$${\mu }_i=P(X_0=s_i),{\mu }_i\ge 0,\sum _{i\in S}{\mu }_i=1$$

有限维分布

给定初始分布 $\mu$ 和转移核 $\mathbf{P}$，马尔可夫链的所有有限维分布由以下公式给出：

$$P(X_0=i_0,X_1=i_1,\dots ,X_n=i_n)={\mu }_{i_0}\cdot p_{i_0,i_1}\cdot p_{i_1,i_2}\cdots p_{i_{n-1},i_n}$$

这一性质称为马尔可夫链的链法则，它表明联合分布完全由初始分布和转移概率决定。

---

n步转移概率

Chapman-Kolmogorov方程

从状态 $i$ 经过 $n$ 步到达状态 $j$ 的概率记为 $p_{ij}^{(n)}$。根据Chapman-Kolmogorov方程：

$$p_{ij}^{(n+m)}=\sum _{k\in S}{p}_{ik}^{(n)}\cdot p_{kj}^{(m)}$$

用矩阵形式表示：

$${\mathbf{P}}^{(n+m)}={\mathbf{P}}^{(n)}\cdot {\mathbf{P}}^{(m)}$$

特别地，$n$ 步转移概率矩阵：

$${\mathbf{P}}^{(n)}={\mathbf{P}}^n$$

计算示例（Python）

import numpy as np
 
# 天气模型的转移矩阵
P = np.array([[0.7, 0.3],
              [0.4, 0.6]])
 
# 计算10步后的转移概率
P_10 = np.linalg.matrix_power(P, 10)
print("10步转移矩阵:")
print(P_10)
# 输出将接近平稳分布的行
 
# 验证Chapman-Kolmogorov方程: P^5 = P^2 @ P^3
P_2 = np.linalg.matrix_power(P, 2)
P_3 = np.linalg.matrix_power(P, 3)
P_5 = np.linalg.matrix_power(P, 5)
print("\n验证 C-K 方程 (P^2 @ P^3 - P^5):")
print(P_2 @ P_3 - P_5)  # 接近零矩阵

---

马尔可夫链的收敛理论

不可约性

一个马尔可夫链称为不可约的（Irreducible），如果对任意两个状态 $i,j\in S$，存在 $n\ge 0$ 使得 $p_{ij}^{(n)}>0$。这意味着：

$$\forall i,j\in S:i\leftrightarrow j$$

不可约性保证了链可以从任意状态到达任意其他状态。

周期性

状态 $i$ 的周期（Period）$d_i$ 定义为：

$$d_i=\gcd \{n\ge 1:p_{ii}^{(n)}>0\}$$

• 若 $d_i=1$，状态是非周期的（Aperiodic）
• 若 $d_i>1$，状态是周期的（Periodic）

周期性质：若马尔可夫链不可约，则所有状态具有相同的周期。

遍历性

状态 $i$ 称为：

• 常返的（Recurrent），如果从 $i$ 出发返回 $i$ 的概率为1：

  $$P_i({\tau }_i^{+}<\infty )=1$$

  其中 ${\tau }_i^{+}=\inf \{n\ge 1:X_n=i\}$ 是首次返回时间。

• 正常返的（Positive Recurrent），如果它是常返的且平均返回时间有限：

  $$E_i[{\tau }_i^{+}]<\infty$$

• 零常返的（Null Recurrent），如果它是常返的但平均返回时间无限。

• 瞬态的（Transient），如果从 $i$ 出发存在正概率永不返回 $i$。

平稳分布的存在性与唯一性

定理（平稳分布存在性）：每个不可约、非周期、正常返的有限马尔可夫链都有唯一的平稳分布。

对于一般状态空间，刘维尔-弗罗贝尼乌斯定理给出：

定理（不动点定理）：设 $\mathbf{P}$ 是不可约转移矩阵，则存在唯一的概率向量 $\pi$ 满足：

$$\pi \mathbf{P}=\pi ,\sum _{i\in S}{\pi }_i=1$$

$\pi$ 称为平稳分布（Stationary Distribution），其分量：

$${\pi }_j=\underset{n\to \infty }{\lim }{p}_{ij}^{(n)}$$

对所有 $i\in S$ 都相同（与初始状态无关）。

细致平衡方程

若存在分布 $\pi$ 满足细致平衡条件（Detailed Balance Equation）：

$${\pi }_i{p}_{ij}={\pi }_j{p}_{ji},\forall i,j\in S$$

则该链是可逆的，且 $\pi$ 是其平稳分布。

收敛速率与谱分析

马尔可夫链的收敛速率与转移矩阵的谱性质密切相关。

谱分解：设 $\mathbf{P}$ 是有限状态空间马尔可夫链的转移矩阵，则 $\mathbf{P}$ 可以分解为：

$$\mathbf{P}=\mathbf{Q\Lambda }{\mathbf{Q}}^{-1}$$

其中 $\mathbf{\Lambda }=\text{diag}({\lambda }_1,{\lambda }_2,\dots ,{\lambda }_n)$ 是特征值对角矩阵。

谱隙（Spectral Gap）：

$${\gamma }^{\ast }=1-{\lambda }_2$$

其中 ${\lambda }_2$ 是第二大特征值（按绝对值）。

收敛界：对于不可约、非周期、正常返的马尔可夫链：

$$\Vert {\mathbf{P}}^n(x,\cdot )-\pi {\Vert }_{\text{TV}}\le C\cdot ({\lambda }^{\ast }{)}^n$$

其中 $\Vert \cdot {\Vert }_{\text{TV}}$ 是全变差距离，${\lambda }^{\ast }=\max (|{\lambda }_2|,|{\lambda }_n|)$ 是第二大特征值的绝对值。

import numpy as np
 
# 分析转移矩阵的谱性质
P = np.array([[0.7, 0.3],
              [0.4, 0.6]])
 
eigenvalues, eigenvectors = np.linalg.eig(P.T)  # 注意转置
eigenvalues = np.real(eigenvalues)
idx = np.argsort(eigenvalues)[::-1]
eigenvalues = eigenvalues[idx]
 
print("特征值:", eigenvalues)
print("谱隙 (1 - |λ₂|):", 1 - abs(eigenvalues[1]))
 
# 验证平稳分布
pi = np.array([4/7, 3/7])
print("平稳分布:", pi)
print("πP =", pi @ P)

---

连续时间马尔可夫链

基本定义

连续时间马尔可夫链（Continuous-Time Markov Chain, CTMC）$\{X(t){\}}_{t\ge 0}$ 满足：

$$P(X(t_{n+1})=j\mid X(t_n)=i,X(t_{n-1}),\dots )=P(X(t_{n+1})=j\mid X(t_n)=i)$$

连续时间链的特征是：在状态 $i$ 停留的时间是指数分布的随机变量。

Kolmogorov微分方程

设 $p_{ij}(t)=P(X(t)=j\mid X(0)=i)$，则满足Kolmogorov向前方程：

$$p_{ij}^{\prime }(t)=\sum _{k\ne j}{p}_{ik}(t)q_{kj}-p_{ij}(t)\sum _{k\ne j}{q}_{jk}$$

矩阵形式：

$${\mathbf{P}}^{\prime }(t)=\mathbf{P}(t)\mathbf{Q}=\mathbf{QP}(t)$$

其中 $\mathbf{Q}$ 是生成元（Infinitesimal Generator）：

$$q_{ij}=\{\begin{array}{ll}{\lambda }_i\cdot r_{ij}& i\ne j\\ -{\lambda }_i& i=j\end{array}$$

这里 ${\lambda }_i$ 是离开状态 $i$ 的总率，$r_{ij}$ 是从 $i$ 转移到 $j$ 的条件概率。

出生-死亡过程

出生-死亡过程（Birth-Death Process）是一种特殊的连续时间马尔可夫链，状态空间为 ${\mathbb{Z}}_{+}$，转移只发生在相邻状态之间。

转移率：

• 出生率：${\lambda }_n$（从状态 $n$ 到 $n+1$）
• 死亡率：${\mu }_n$（从状态 $n$ 到 $n-1$）

Kolmogorov平衡方程：

$${\lambda }_{n-1}{\pi }_{n-1}+{\mu }_{n+1}{\pi }_{n+1}=({\lambda }_n+{\mu }_n){\pi }_n$$

平稳分布（通过递归可得）：

$${\pi }_n={\pi }_0\prod _{k=0}^{n-1}\frac{{\lambda }_k}{{\mu }_{k+1}},n\ge 1$$

Poisson过程

Poisson过程是连续时间马尔可夫链的经典例子。

定义：计数过程 $\{N(t){\}}_{t\ge 0}$ 称为参数为 $\lambda$ 的Poisson过程，如果：

1. $N(0)=0$
2. 独立增量：$N(t)-N(s)$ 与 $N(v)-N(u)$ 独立（$s>v\ge u>t$）
3. 平稳增量：$N(t)-N(s)\sim \text{Poisson}(\lambda (t-s))$
4. 轨道右连续且为阶梯函数

转移率矩阵：

$$q_{ij}=\{\begin{array}{ll}\lambda & j=i+1\\ -\lambda & j=i\\ 0& \text{其他}\end{array}$$

import numpy as np
import matplotlib.pyplot as plt
 
def poisson_process(lam, T):
    """模拟参数为λ的Poisson过程"""
    times = [0]
    counts = [0]
    t = 0
    n = 0
    
    while t < T:
        # 指数分布的等待时间
        wait = np.random.exponential(1 / lam)
        t += wait
        if t > T:
            break
        n += 1
        times.append(t)
        counts.append(n)
    
    # 阶梯函数需要水平部分
    times_full = times + [T]
    counts_full = counts + [n]
    
    return times_full, counts_full
 
# 模拟并绘图
lam = 5
T = 2
times, counts = poisson_process(lam, T)
 
plt.figure(figsize=(10, 4))
plt.step(times, counts, where='post', linewidth=2)
plt.xlabel('时间 t')
plt.ylabel('计数 N(t)')
plt.title(f'Poisson过程 (λ = {lam})')
plt.grid(True, alpha=0.3)
plt.savefig('/tmp/poisson_process.png', dpi=150)

---

稳态分布

平稳分布的定义

若一个概率分布 $\pi$ 满足：

$$\pi \mathbf{P}=\pi$$

则称 $\pi$ 为该马尔可夫链的平稳分布（Stationary Distribution）。

极限分布

对于许多马尔可夫链，当 $n\to \infty$ 时，$n$ 步转移概率 $P(X_n=j)$ 趋向于一个与初始状态无关的极限分布 ${\pi }_j$。

遍历定理

遍历定理（Ergodic Theorem）指出：对于不可约、非周期、正常返的马尔可夫链，时间平均趋向于空间平均。

即对任意函数 $f:S\to \mathbb{R}$：

$$\frac{1}{n}\sum _{t=0}^{n-1}f(X_t)\stackrel{n\to \infty }{\to }\sum _{i\in S}{\pi }_{i}f(i)\text{（几乎必然）}$$

---

状态分类

可达与互通

• 可达：若存在 $n\ge 0$ 使得 $p_{ij}^{(n)}>0$，则称 $j$ 可从 $i$ 可达，记为 $i\to j$
• 互通：若 $i\to j$ 且 $j\to i$，则称 $i$ 与 $j$ 互通，记为 $i\leftrightarrow j$

周期性

状态 $i$ 的周期 $d$ 定义为：

$$d=\gcd \{n\ge 1:p_{ii}^{(n)}>0\}$$

• 若 $d=1$，则状态是非周期的
• 若 $d>1$，则状态是周期的

常返与瞬态

• 常返（Recurrent）：从状态 $i$ 出发，以概率1返回 $i$
• 瞬态（Transient）：从状态 $i$ 出发，存在正概率永不返回 $i$

---

马尔可夫链蒙特卡洛方法（MCMC）

马尔可夫链的一个强大应用是马尔可夫链蒙特卡洛方法（Markov Chain Monte Carlo），用于从复杂分布中采样。

Metropolis-Hastings算法

Metropolis-Hastings算法的核心思想是构造一个以目标分布 $\pi$ 为平稳分布的马尔可夫链。

算法步骤：

1. 初始化 $X_0$，设置提议分布 $q(x^{\prime }\mid x)$
2. 对 $t=0,1,2,\dots$：
   • 从 $q(\cdot \mid X_t)$ 采样候选状态 $x^{\prime }$
   • 计算接受率： $$\alpha =\min \left(1,\frac{\pi (x^{\prime })q(X_t\mid x^{\prime })}{\pi (X_t)q(x^{\prime }\mid X_t)}\right)$$
   • 以概率 $\alpha$ 接受：$X_{t+1}=x^{\prime }$，否则拒绝：$X_{t+1}=X_t$

import numpy as np
 
def metropolis_hastings_mcmc(pdf, proposal_std, n_samples, dim=2):
    """
    Metropolis-Hastings算法采样
    
    Args:
        pdf: 目标分布（未归一化）
        proposal_std: 提议分布的标准差
        n_samples: 采样数量
        dim: 维度
    
    Returns:
        samples: 采样序列
    """
    samples = np.zeros((n_samples, dim))
    current = np.random.randn(dim) * 0.5  # 初始点
    
    for i in range(n_samples):
        # 提议分布：各向同性正态分布
        proposal = current + np.random.randn(dim) * proposal_std
        
        # 计算接受率
        alpha = min(1, pdf(proposal) / pdf(current))
        
        # 接受/拒绝
        if np.random.random() < alpha:
            current = proposal
        
        samples[i] = current
    
    return samples
 
# 示例：采样二元正态分布
def bivariate_normal_pdf(x):
    mean = [0, 0]
    cov = [[1, 0.6], [0.6, 1]]
    det = np.linalg.det(cov)
    inv_cov = np.linalg.inv(cov)
    diff = x - mean
    return np.exp(-0.5 * diff @ inv_cov @ diff) / (2 * np.pi * np.sqrt(det))
 
samples = metropolis_hastings_mcmc(bivariate_normal_pdf, 0.5, 10000)
print(f"采样均值: {samples.mean(axis=0)}")
print(f"采样协方差:\n{np.cov(samples.T)}")

Gibbs采样

Gibbs采样是Metropolis-Hastings的特殊情况，当提议分布为条件分布时，接受率为1。

全条件分布：给定其他变量时单个变量的条件分布：

$$X_i^{(t+1)}\sim p\left(x_i\mid x_1^{(t+1)},\dots ,x_{i-1}^{(t+1)},x_{i+1}^{(t)},\dots ,x_d^{(t)}\right)$$

def gibbs_sampling_gaussian(mean, cov, n_samples):
    """
    对多元高斯分布进行Gibbs采样
    使用坐标交替更新
    """
    d = len(mean)
    samples = np.zeros((n_samples, d))
    x = np.zeros(d)  # 初始值
    
    for i in range(n_samples):
        for j in range(d):
            # 计算条件分布的均值和方差
            sigma_j_sq = 1.0 / cov[j, j]
            mu_j = mean[j] + sigma_j_sq * sum(
                cov[j, k] * (x[k] - mean[k]) for k in range(d) if k != j
            )
            x[j] = np.random.normal(mu_j, np.sqrt(sigma_j_sq))
        samples[i] = x
    
    return samples

收敛诊断

MCMC链的收敛诊断是一个重要但困难的问题。

常用方法：

1. Trace Plot（迹图）：观察链是否”平稳”
2. Geweke诊断：比较链前后部分的均值是否一致
3. Gelman-Rubin统计量（$\hat{R}$）：比较多个链的方差

def gelman_rubin(chains):
    """
    Gelman-Rubin诊断
    
    Args:
        chains: list of arrays, 每个chain是独立的马尔可夫链
    
    Returns:
        R_hat: 潜在尺度缩减因子，接近1表示收敛
    """
    m = len(chains)  # 链的数量
    n = chains[0].shape[0]  # 每个链的样本数
    
    # 计算每个链的均值和方差
    chain_means = np.array([np.mean(chain, axis=0) for chain in chains])
    chain_vars = np.array([np.var(chain, axis=0, ddof=1) for chain in chains])
    
    # 链间方差
    B = n / (m - 1) * np.sum((chain_means - np.mean(chain_means, axis=0))**2, axis=0)
    
    # 链内方差
    W = np.mean(chain_vars, axis=0)
    
    # 总体方差估计
    var_hat = (n - 1) / n * W + B / n
    
    # R_hat
    R_hat = np.sqrt(var_hat / W)
    
    return R_hat
 
def trace_plot(samples, name="Chain"):
    """绘制迹图"""
    plt.figure(figsize=(12, 4))
    plt.plot(samples, alpha=0.7)
    plt.axhline(np.mean(samples), color='red', linestyle='--', label='Mean')
    plt.xlabel('Iteration')
    plt.ylabel('Value')
    plt.title(f'{name} Trace Plot')
    plt.legend()
    plt.grid(True, alpha=0.3)

收敛标准：

• $\hat{R}<1.1$：通常认为链已收敛
• $\hat{R}<1.01$：更严格的标准

---

应用实例：PageRank

Google的PageRank算法本质上就是一个马尔可夫链。

基本思想

• 网页的重要性取决于指向它的网页数量和质量
• 用户随机浏览网页的行为可以用马尔可夫链建模
• PageRank值就是该马尔可夫链的平稳分布

数学表述

设网页 $i$ 的PageRank为 $PR(i)$，$L_j$ 为页面 $j$ 的出链数，则：

$$PR(i)=\frac{1-d}{N}+d\sum _{j\in B_i}\frac{PR(j)}{L_j}$$

其中 $d$ 是阻尼因子（通常取0.85），$B_i$ 是指向 $i$ 的页面集合。

PageRank与马尔可夫链

PageRank等价于随机冲浪者（Random Surfer）模型的平稳分布：

1. 以概率 $1-d$ 随机跳转到任意页面
2. 以概率 $d$ 按链接随机跳转

def pagerank(adjacency_matrix, d=0.85, tol=1e-6):
    """
    计算PageRank
    
    Args:
        adjacency_matrix: 邻接矩阵 (N x N)
        d: 阻尼因子
        tol: 收敛阈值
    
    Returns:
        pagerank_scores: PageRank分数
    """
    n = adjacency_matrix.shape[0]
    
    # 归一化为转移矩阵
    out_degree = adjacency_matrix.sum(axis=1, keepdims=True)
    out_degree[out_degree == 0] = 1  # 处理悬挂节点
    P = adjacency_matrix / out_degree
    
    # 初始分布
    pi = np.ones(n) / n
    
    # 幂迭代
    while True:
        pi_new = (1 - d) / n + d * (pi @ P)
        if np.linalg.norm(pi_new - pi, 1) < tol:
            break
        pi = pi_new
    
    return pi

---

应用实例：天气模型

考虑一个简化的天气模型：

状态	晴	雨
晴	0.7	0.3
雨	0.4	0.6

转移概率矩阵：

$$\mathbf{P}=\left(\begin{array}{cc}0.7& 0.3\\ 0.4& 0.6\end{array}\right)$$

稳态分布计算

设稳态分布为 $({\pi }_1,{\pi }_2)$，满足：

$$\{\begin{array}{l}{\pi }_1=0.7{\pi }_1+0.4{\pi }_2\\ {\pi }_1+{\pi }_2=1\end{array}$$

解得：${\pi }_1=\frac{4}{7}\approx 0.571$，${\pi }_2=\frac{3}{7}\approx 0.429$

这意味着长期来看，约57.1%的时间是晴天，42.9%的时间是雨天。

模拟验证

import numpy as np
 
def simulate_weather(n_steps, transition_matrix, initial_state=0):
    """
    模拟天气马尔可夫链
    
    Args:
        n_steps: 模拟步数
        transition_matrix: 转移概率矩阵
        initial_state: 初始状态 (0=晴, 1=雨)
    
    Returns:
        states: 状态序列
    """
    states = [initial_state]
    current = initial_state
    
    for _ in range(n_steps - 1):
        probs = transition_matrix[current]
        current = np.random.choice(len(probs), p=probs)
        states.append(current)
    
    return np.array(states)
 
# 参数
P = np.array([[0.7, 0.3],
              [0.4, 0.6]])
 
# 模拟10000步
np.random.seed(42)
states = simulate_weather(10000, P)
 
# 计算经验频率
sunny_freq = np.mean(states == 0)
rainy_freq = np.mean(states == 1)
 
print(f"晴天频率: {sunny_freq:.4f}")
print(f"雨天频率: {rainy_freq:.4f}")
print(f"理论晴天概率: {4/7:.4f}")
print(f"理论雨天概率: {3/7:.4f}")

---

与其他模型的关系

隐马尔可夫模型（HMM）

隐马尔可夫模型是马尔可夫链的扩展，其中状态是”隐藏”的，只能通过观测序列来推断。详见 隐马尔可夫模型。

条件随机场（CRF）

条件随机场是一种判别式模型，用于序列标注任务。相比HMM，CRF可以建模更复杂的依赖关系。详见后续的概率图模型文档。

强化学习中的马尔可夫决策过程

马尔可夫决策过程（MDP）是马尔可夫链的决策扩展，加入了动作和奖励。详见 MDP。

---

参考

Footnotes

1. 本文档参考了《概率论与数理统计》相关章节和斯坦福CS228课程材料 ↩