2.8 随机过程与概率图模型

前面几节处理的对象是有限个随机变量。2.7 节讨论的极限定理研究一列独立同分布随机变量的累积行为：大数定律说明样本均值依概率收敛到期望，中心极限定理说明标准化和依分布收敛到标准正态分布。这些结论都以指标之间相互独立为前提，样本之间不存在时序或结构上的依赖。实际数据很少满足这一点：今天的股价与昨天的收盘价相关，句子中的每个词依赖它前面的词，社交网络中一个用户的活跃度受其好友的影响。这些对象由一列或一组相互依赖的随机变量描述，需要一套新的语言。

随机过程把随机变量按某个指标集排列，用联合分布刻画指标之间的依赖结构。指标集可以表示时间，也可以表示空间位置或图上的节点。依赖结构一旦固定下来，原本需要逐一写出全部联合分布的问题，就转化为少数几个参数或图结构的问题。本节先给出随机过程的定义与基本分类，再依次讨论状态离散的马尔可夫链、隐马尔可夫模型与泊松过程，取值连续的布朗运动与高斯过程，最后把变量之间的依赖关系用图表示，引入概率图模型与因果图。

这一章的内容是统计推断与序列建模的共同基础。马尔可夫链是 MCMC 采样的引擎，隐马尔可夫模型是序列标注与语音识别的经典工具，高斯过程给出贝叶斯非参数的回归框架，概率图模型统一了贝叶斯网络与马尔可夫随机场的因子分解。读懂这些模型，才能在后续章节中处理带结构的高维数据。

2.8.1 随机过程的定义与分类

随机变量族与指标集

随机过程（stochastic process）是定义在概率空间 $(\mathrm{\Omega },\mathcal{F},P)$ 上的一族随机变量 $\{X_t{\}}_{t\in T}$，其中 $T$ 称为指标集（index set），每个 $t\in T$ 对应一个随机变量 $X_t$。固定 $\omega \in \mathrm{\Omega }$ 时，映射 $t\mapsto X_t(\omega )$ 是定义在 $T$ 上的一个函数，称为该过程的一条样本轨道（sample path）或一次实现。同一过程在每个 $\omega$ 上给出一条轨道，全部轨道构成过程的全部不确定性。

指标集 $T$ 通常取 $\text{N}$、$\text{Z}$、$[0,+\mathrm{\infty })$ 或 $\text{R}$。$T$ 离散时称为离散时间过程，习惯上写成 $X_1,X_2,\dots$ 或 $\{X_n{\}}_{n\ge 0}$；$T$ 连续时称为连续时间过程，写成 $\{X_t{\}}_{t\ge 0}$。每个 $X_t$ 的取值集合称为状态空间（state space），记作 $S$。状态空间可以是有限集、可数集或 ${\text{R}}^d$。指标集与状态空间的组合给出了四类常见情形，它们的处理工具差别很大。

指标集 $T$	状态空间 $S$	典型过程
离散	可数	随机游走、马尔可夫链
离散	连续	时间序列、语音信号
连续	可数	泊松过程、计数过程
连续	连续	布朗运动、高斯过程

有限维分布族

随机过程在数学上由它的全部有限维分布决定。对任意 $n\ge 1$ 与任意时刻 $t_1<t_2<\cdots <t_n$，联合分布

$$F_{t_1,\dots ,t_n}(x_1,\dots ,x_n)=P(X_{t_1}\le x_1,\dots ,X_{t_n}\le x_n)$$

称为过程的一个 $n$ 维边缘分布。全部这类联合分布组成的集合称为有限维分布族（family of finite-dimensional distributions）。两个随机过程若有限维分布族相同，则在统计上无法区分，因此建模时只需给出这个分布族。

反过来，给定一族满足相容性条件的有限维分布，可以构造出概率空间与随机过程以其为有限维分布。这一结论称为柯尔莫哥洛夫扩张定理（Kolmogorov extension theorem）。相容性要求低维分布是高维分布取边缘的结果，同时分布对时刻的排序保持不变。扩张定理把描述一个过程的问题转化为描述一族相容的联合分布，也为连续时间过程的存在性提供了保证。

矩函数

完整描述有限维分布族往往不现实，实践中常用前两阶矩来概括过程的结构。均值函数定义为 $m(t)=E[X_t]$，它给出过程在每个时刻的平均水平。自协方差函数定义为

$$C(s,t)=\mathrm{Cov}(X_s,X_t)=E[(X_s-m(s))(X_t-m(t))]$$

它刻画不同时刻取值之间的线性相关程度。当 $s=t$ 时，$C(t,t)=\mathrm{Var}(X_t)$ 退化为方差。均值函数与自协方差函数合称二阶矩结构，高斯过程的全部有限维分布由它们完全确定，这是高斯过程便于处理的原因之一。

典型例子

随机游走（random walk）是最简单的依赖序列。设 $X_1,X_2,\dots$ 独立同分布，$P(X_i=1)=p$，$P(X_i=-1)=1-p$，部分和 $S_0=0$，$S_n=X_1+\cdots +X_n$ 构成随机游走。状态空间为 $\text{Z}$，指标集为 $\text{N}$。$p=1/2$ 时称为简单对称随机游走，对任意 $s\le t$ 有 $\mathrm{Cov}(S_s,S_t)=s$，两个时刻的相关性只依赖较早的时刻。随机游走是布朗运动的离散版本，也是马尔可夫链的原型。

计数过程（counting process）记录到时刻 $t$ 为止发生的事件总数，记作 $N_t$。这类过程满足 $N_0=0$，取值非负整数，轨道右连续且单调不减，每发生一次事件就跳 1。排队系统到达的顾客数、放射性衰变中探测器记录到的粒子数、网站在单位时间内的访问次数都属于计数过程。泊松过程是计数过程中最重要的一族。

信号序列是取值连续的时间序列，如语音波形、气温序列、股票价格序列。这类过程的状态空间为 $\text{R}$ 或 ${\text{R}}^d$，相邻时刻的取值通过某种平滑或相关性结构联系起来。平稳过程与谱分析为这类数据提供了标准工具，在 2.8.8 节讨论。

有限维分布族的作用

定义随机过程时容易只关注单个时刻的分布，忽略不同时刻之间的联合结构。判断两个过程是否相同、过程是否平稳、增量是否独立，都要用到有限维分布。扩张定理说明有限维分布族与过程本身等价，这为后续所有基于联合分布的定义提供了依据。实际建模时很少有人写出全部有限维分布，通常只指定均值、协方差或转移概率，由它们隐含地确定分布族。

import numpy as np

rng = np.random.default_rng(0)

# 简单对称随机游走的一条样本轨道
n_steps = 500
steps = rng.choice([-1, 1], size=n_steps)
path = np.concatenate([[0], np.cumsum(steps)])
print(f"轨道长度 = {path.shape[0]}")   # 501

# 用大量轨道验证 Var(S_n) = n
n_paths = 20000
final = rng.choice([-1, 1], size=(n_paths, n_steps)).sum(axis=1)
print(f"模拟方差 = {final.var():.2f}，理论值 n = {n_steps}")
print(f"模拟标准差 = {final.std():.2f}，理论值 sqrt(n) = {np.sqrt(n_steps):.2f}")
# 模拟方差 = 498.93，理论值 n = 500
# 模拟标准差 = 22.34，理论值 sqrt(n) = 22.36

2.8.2 马尔可夫链

马尔可夫性质

马尔可夫链（Markov chain）是状态空间可数、指标集离散的随机过程，它的依赖结构由一条简洁的假设规定。设 $\{X_n{\}}_{n\ge 0}$ 取值于可数集 $S$，若对任意 $n$ 与任意状态 $i,j,i_0,\dots ,i_{n-1}$，只要条件概率有定义，就有

$$P(X_{n+1}=j\mid X_n=i,X_{n-1}=i_{n-1},\dots ,X_0=i_0)=P(X_{n+1}=j\mid X_n=i)$$

则称该过程具有马尔可夫性质（Markov property），也称无后效性。条件的含义是：给定当前状态后，未来的分布与更早的历史无关，全部历史信息都被压缩进了当前状态。在马尔可夫性质下，条件概率 $P(X_{n+1}=j\mid X_n=i)$ 称为转移概率（transition probability），记作 $p_{ij}(n)$。当它与 $n$ 无关时，称链是时齐的，记作 $p_{ij}$，本节默认讨论时齐链。

转移矩阵与初始分布

把全部转移概率排成矩阵 $P=(p_{ij}{)}_{i,j\in S}$，称为转移矩阵（transition matrix）。它是随机矩阵：每个元素非负，每行之和为 1。随机矩阵的第 $i$ 行给出从状态 $i$ 出发下一步的分布。初始时刻的分布 ${\pi }_0=({\pi }_0(i){)}_{i\in S}$ 给出 ${\pi }_0(i)=P(X_0=i)$。初始分布与转移矩阵合起来确定了整条链的全部有限维分布。一行 $n$ 步后的状态分布由 ${\pi }_n={\pi }_0{P}^n$ 给出，$P^n$ 的第 $(i,j)$ 项 $p_{ij}^{(n)}$ 表示从 $i$ 出发经过 $n$ 步到达 $j$ 的概率。

$n$ 步转移与查普曼-科尔莫戈罗夫方程

多步转移概率满足查普曼-科尔莫戈罗夫方程（Chapman-Kolmogorov equation）：对任意 $m,n\ge 0$，

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

方程的含义是把 $m+n$ 步的转移按中间时刻的状态 $k$ 展开，先走 $m$ 步到 $k$，再走 $n$ 步到 $j$，对所有中间状态求和。用矩阵形式书写，方程即 $P^{m+n}=P^m{P}^n$，与矩阵乘法的定义一致。方程说明多步转移完全由一步转移决定，链的长期行为因此可以从单个矩阵 $P$ 读出。

状态分类

链的长期行为取决于各个状态的返回性质。对状态 $i$，令 $f_i$ 为从 $i$ 出发最终返回 $i$ 的概率。若 $f_i=1$，称 $i$ 为常返态（recurrent state）；若 $f_i<1$，称 $i$ 为瞬态（transient state）。常返态又按时均返回时间分为正常返与零常返，正常返指平均返回时间有限。有限状态的链中不存在零常返态。若两个状态可以相互到达，称它们互通，互通关系把状态空间划分为若干互通类。

状态的周期（period）定义为使 $p_{ii}^{(n)}>0$ 的步数 $n$ 的最大公约数，记作 $d(i)$。$d(i)=1$ 时称该状态非周期。周期意味着返回只能发生在某些固定的步数上，例如两状态链在两状态间来回跳动，每个状态的周期为 2。若某状态满足 $p_{ii}=1$，称它为吸收态（absorbing state），一旦进入就不再离开。吸收态常返且周期为 1。若链只有一个互通类且非周期，称链不可约且非周期。

平稳分布与细致平衡

分布 $\pi =(\pi (i){)}_{i\in S}$ 若满足 $\pi \ge 0$、$\sum _i\pi (i)=1$ 以及

$$\pi =\pi P,\text{即}\pi (j)=\sum _{i\in S}\pi (i)p_{ij}$$

则称 $\pi$ 为链的平稳分布（stationary distribution）。若初始分布取 $\pi$，则任意时刻的边际分布都是 $\pi$，链在分布意义下保持不变。平稳分布是 $P$ 的左特征向量，对应特征值 1 的特征空间。

一个更强的条件是细致平衡（detailed balance）：对任意 $i,j$，

$$\pi (i)p_{ij}=\pi (j)p_{ji}$$

满足细致平衡的链称为可逆链。把细致平衡两边对 $i$ 求和，左边得到 $\sum _i\pi (i)p_{ij}=(\pi P{)}_j$，右边得到 $\pi (j)\sum _i{p}_{ji}=\pi (j)$，因此细致平衡蕴含平稳性。细致平衡的条件比平稳性强，验证起来却更方便，MCMC 算法正是通过设计满足细致平衡的转移规则来构造以目标分布为平稳分布的链。

遍历性与极限分布

当链不可约、非周期且正常返时，存在唯一的平稳分布 $\pi$，且对任意初始状态 $i$，

$$p_{ij}^{(n)}⟶\pi (j),n\to \mathrm{\infty }$$

这一性质称为遍历性（ergodicity），此时的平稳分布也称为极限分布。遍历性的含义是长期分布与初始条件无关，链会忘记起点。当链是周期链时，$p_{ij}^{(n)}$ 不收敛，但沿时间平均仍趋于平稳分布。遍历定理进一步保证，对任意函数 $f$，时间平均依概率收敛到空间平均：

$$\frac{1}{n}\sum _{k=0}^{n-1}f(X_k)⟶\sum _{i\in S}f(i)\pi (i)$$

这一定理是 MCMC 方法的理论依据，它允许用一条长轨道的时间平均来近似目标分布的期望。

平稳分布与极限分布的区别

平稳分布只要求初始分布与转移后分布相同，周期链也可以有平稳分布。极限分布要求 $p_{ij}^{(n)}$ 逐项收敛，需要非周期条件。二者在不可约、非周期、正常返的链上重合。设计采样算法时关心的是平稳分布，用轨道平均替代空间平均时依赖遍历定理，两者都需要链的不可约与常返性质。

应用接口：PageRank 与 MCMC

PageRank 把网页排序问题转化为马尔可夫链的平稳分布问题。把网页看作状态，超链接看作转移，转移概率按出链数均分。若某网页没有出链，规定它等概率跳到任意网页；为避免链被少数网页的链接结构困住，再以概率 $1-d$ 跳到全网页上的均匀分布，得到

$$P_{ij}=d\cdot \frac{1}{\text{out}(i)}+\frac{1-d}{N},d\approx 0.85$$

修正后的链不可约、非周期，平稳分布存在且唯一，第 $i$ 个分量就是页面 $i$ 的排名权重。实际计算用幂迭代 ${\pi }_{k+1}={\pi }_{k}P$ 反复右乘转移矩阵，直到相邻两次结果接近，这与求解左特征向量的幂法一致。

马尔可夫链蒙特卡洛（Markov chain Monte Carlo，MCMC）用马尔可夫链生成目标分布的样本。目标是构造一条以目标分布 $\pi$ 为平稳分布的链，然后沿轨道采样。Metropolis-Hastings 算法从当前状态按提议分布生成候选状态，以接受概率

$$\alpha (i,j)=min\{1,\frac{\pi (j)q(j,i)}{\pi (i)q(i,j)}\}$$

决定是否接受，其中 $q$ 是提议分布。接受概率的设计使链满足细致平衡，$\pi$ 成为平稳分布。Gibbs 采样是 Metropolis-Hastings 的特例，它轮流从各变量的条件分布采样，接受概率恒为 1。MCMC 处理的是无法直接采样的高维分布，在 2.12 节讨论贝叶斯计算时会详细展开。

import numpy as np

# 转移矩阵：三状态天气模型（0 晴，1 阴，2 雨）
P = np.array([[0.7, 0.2, 0.1],
              [0.3, 0.4, 0.3],
              [0.2, 0.3, 0.5]])

# 幂迭代求解平稳分布：pi_{k+1} = pi_k P
pi = np.array([1 / 3, 1 / 3, 1 / 3])
for _ in range(200):
    pi = pi @ P
print(f"平稳分布（幂迭代）= {np.round(pi, 6)}")
# [0.432432 0.270270 0.297297]

# 交叉验证：P^T 对应特征值 1 的左特征向量
vals, vecs = np.linalg.eig(P.T)
idx = np.argmin(np.abs(vals - 1.0))
pi_eig = np.real(vecs[:, idx])
pi_eig = pi_eig / pi_eig.sum()
print(f"平稳分布（特征向量）= {np.round(pi_eig, 6)}")

# 模拟验证：从状态 0 出发跑长轨道，统计各状态频率
rng = np.random.default_rng(7)
state = 0
counts = np.zeros(3)
n_steps = 200_000
for _ in range(n_steps):
    state = rng.choice(3, p=P[state])
    counts[state] += 1
print(f"模拟频率 = {np.round(counts / n_steps, 4)}")
# [0.4325 0.2709 0.2966]，与平稳分布一致

2.8.3 隐马尔可夫模型

状态与观测的双层结构

隐马尔可夫模型（hidden Markov model，HMM）把观测序列背后的状态序列建模为马尔可夫链，观测由状态按发射分布生成。模型包含两部分：状态序列 $\{X_t\}$ 是不可观测的马尔可夫链，转移概率为 $a_{ij}=P(X_{t+1}=j\mid X_t=i)$；观测序列 $\{Y_t\}$ 在每个时刻 $t$ 由当前状态 $X_t$ 按发射概率 $b_i(y)=P(Y_t=y\mid X_t=i)$ 生成。给定状态序列时，各观测相互独立。模型的参数由初始分布 $\pi$、转移矩阵 $A=(a_{ij})$ 与发射矩阵 $B=(b_i(y))$ 组成，记作三元组 $\lambda =(\pi ,A,B)$。

双层结构对应许多实际场景。在语音识别中，隐状态表示音素，观测是声学特征；在序列标注中，隐状态表示词性，观测是词；在手势识别中，隐状态表示动作类型，观测是传感器读数。要使用模型，需要解决三个问题：给定参数与观测序列，计算该序列出现的概率（评估问题）；给定参数与观测序列，推断最可能的状态序列（解码问题）；给定观测序列，估计参数（学习问题）。前向算法、维特比算法与鲍姆-韦尔奇算法分别解决这三个问题。

前向算法

前向算法（forward algorithm）计算观测序列的概率 $P(Y_1,\dots ,Y_T)$。定义前向变量

$${\alpha }_t(i)=P(Y_1=y_1,\dots ,Y_t=y_t,X_t=i\mid \lambda )$$

它表示到时刻 $t$ 为止观测序列恰好为给定值且当前状态为 $i$ 的联合概率。初始化 ${\alpha }_1(i)=\pi (i)b_i(y_1)$，递推关系为

$${\alpha }_{t+1}(j)=[\sum _i{\alpha }_t(i)a_{ij}]b_j(y_{t+1})$$

递推把时刻 $t$ 的所有状态按转移概率汇聚到状态 $j$，再乘以发射概率。最终

$$P(Y_1,\dots ,Y_T\mid \lambda )=\sum _i{\alpha }_T(i)$$

直接枚举状态序列需要 $O(|S{|}^T)$ 次计算，前向算法利用马尔可夫结构把复杂度降到 $O(|S{|}^{2}T)$，这是动态规划在概率计算中的一次典型应用。

后向算法

后向变量定义为

$${\beta }_t(i)=P(Y_{t+1}=y_{t+1},\dots ,Y_T=y_T\mid X_t=i,\lambda )$$

它表示在时刻 $t$ 处于状态 $i$ 的条件下，未来观测序列的概率。初始化 ${\beta }_T(i)=1$，递推关系为

$${\beta }_t(i)=\sum _j{a}_{ij}{b}_j(y_{t+1}){\beta }_{t+1}(j)$$

后向算法从时刻 $T$ 向前递推，与前向算法方向相反。两者结合给出在时刻 $t$ 处于状态 $i$ 的后验概率

$${\gamma }_t(i)=P(X_t=i\mid Y_1,\dots ,Y_T)=\frac{{\alpha }_t(i){\beta }_t(i)}{P(Y_1,\dots ,Y_T)}$$

类似的乘积 ${\xi }_t(i,j)\propto {\alpha }_t(i)a_{ij}{b}_j(y_{t+1}){\beta }_{t+1}(j)$ 给出相邻两次状态转移的后验概率。这两个量是鲍姆-韦尔奇算法的基本构件，也是计算平滑分布的标准工具。

维特比算法

维特比算法（Viterbi algorithm）求最可能的状态序列。定义

$${\delta }_t(i)=\underset{x_1,\dots ,x_{t-1}}{max}P(x_1,\dots ,x_{t-1},X_t=i,Y_1,\dots ,Y_t\mid \lambda )$$

它表示以状态 $i$ 结尾的最优路径的概率。初始化 ${\delta }_1(i)=\pi (i)b_i(y_1)$，递推关系为

$${\delta }_t(j)=[\underset{i}{max}{\delta }_{t-1}(i)a_{ij}]b_j(y_t)$$

同时记录使最大值达到的 $i$，保存在回溯指针 ${\psi }_t(j)$ 中。递推结束后取 ${\delta }_T$ 的最大值位置作为终点，再沿回溯指针恢复整条状态序列。维特比算法把前向算法中的求和换成取最大值，复杂度同样为 $O(|S{|}^{2}T)$。它给出的是一条联合概率最大的路径，这一序列不同于对每个时刻取 ${\gamma }_t(i)$ 最大值得到的序列，后者的拼接可能违背转移约束。

鲍姆-韦尔奇算法

参数未知时用鲍姆-韦尔奇算法（Baum-Welch algorithm）估计 $\lambda$。它属于期望最大化（expectation-maximization，EM）方法，在隐状态不可观测的情形下迭代提升观测序列的似然。算法的每次迭代分两步。期望步用当前参数计算后验概率 ${\gamma }_t(i)$ 与 ${\xi }_t(i,j)$；最大化步用这些后验作为权重更新参数：

$${\hat{a}}_{ij}=\frac{\sum _{t=1}^{T-1}{\xi }_t(i,j)}{\sum _{t=1}^{T-1}{\gamma }_t(i)},{\hat{b}}_i(y)=\frac{\sum _{t:y_t=y}{\gamma }_t(i)}{\sum _{t=1}^T{\gamma }_t(i)},\hat{\pi }(i)={\gamma }_1(i)$$

更新公式的含义是把当前后验看成隐状态的软计数，用这些计数估计转移频率与发射频率。可以证明每次迭代后观测序列的似然不减，算法收敛到似然的局部极大值。初始参数的选取影响收敛结果，实践中常多次随机初始化后取似然最大的一组。

在序列标注与语音识别中的位置

HMM 在 20 世纪后半叶成为语音识别与序列标注的标准模型。声学模型把音素作为隐状态，把每帧声学特征作为观测，用前向算法评分，用维特比算法解码出音素序列，再接语言模型得到词序列。在自然语言处理中，HMM 用于词性标注与命名实体识别，隐状态对应词性标签，观测对应词本身。

HMM 的局限来自两个假设：观测在给定状态下相互独立，且状态之间只存在一阶依赖。观测之间的长距离相关与重叠特征难以用发射分布刻画。条件随机场放宽了观测的独立性假设，改用全局归一化直接建模 $P(Y\mid X)$，在序列标注任务上通常优于 HMM。

import numpy as np

# 两状态 HMM：状态 0 健康、状态 1 患病；观测取 3 个符号
A = np.array([[0.7, 0.3],
              [0.4, 0.6]])            # 转移矩阵
B = np.array([[0.9, 0.05, 0.05],      # 发射矩阵 B[i, k] = P(观测 k | 状态 i)
              [0.2, 0.30, 0.50]])
pi = np.array([0.6, 0.4])             # 初始分布
obs = [0, 1, 2, 0]                    # 观测序列

def forward(A, B, pi, obs):
    T = len(obs)
    n = len(pi)
    alpha = np.zeros((T, n))
    alpha[0] = pi * B[:, obs[0]]
    for t in range(1, T):
        alpha[t] = (alpha[t - 1] @ A) * B[:, obs[t]]
    return alpha

alpha = forward(A, B, pi, obs)
print(np.round(alpha, 6))
print(f"观测序列的概率 P(O) = {alpha[-1].sum():.6f}")
# P(O) = 0.011912

2.8.4 泊松过程

齐次泊松过程的定义

泊松过程（Poisson process）是计数过程中最基本的一族。强度为 $\lambda >0$ 的齐次泊松过程 $\{N_t{\}}_{t\ge 0}$ 满足三条性质：$N_0=0$；过程具有独立增量，即任意不重叠时间区间上的增量相互独立；过程具有平稳增量，即增量只依赖区间长度，$N_{t+s}-N_s$ 服从均值为 $\lambda t$ 的泊松分布，

$$P(N_{t+s}-N_s=k)=\frac{(\lambda t{)}^k}{k!}{e}^{-\lambda t},k=0,1,2,\dots$$

三条性质合起来刻画了事件在时间上完全随机地发生这一情形：任意小区的发生概率与区间长度成正比，两个事件在同一瞬间发生的概率为零。$\lambda$ 表示单位时间内的平均事件数，因此 $E[N_t]=\mathrm{Var}(N_t)=\lambda t$，均值与方差相等是泊松分布的特征，也是检验数据是否适合泊松模型的线索。

到达间隔与指数分布

令 $T_1$ 为第一次事件的发生时刻，$T_k$ 为第 $k-1$ 次与第 $k$ 次事件之间的时间间隔。由泊松分布可推出

$$P(T_1>t)=P(N_t=0)=e^{-\lambda t}$$

因此到达间隔 $T_1,T_2,\dots$ 独立同分布于参数为 $\lambda$ 的指数分布，其密度为 $f(t)=\lambda e^{-\lambda t}$，均值为 $1/\lambda$。指数分布具有无记忆性：$P(T>s+t\mid T>s)=P(T>t)$，已知已经等待了 $s$ 时间，后续等待的分布与重新开始相同。这一性质与泊松过程的平稳增量等价。

反过来，若一列独立同分布的指数间隔，其到达时刻构成齐次泊松过程。同一过程的两种描述在建模时各有方便之处：关心事件计数时用泊松分布，关心等待时间时用指数分布。

非齐次泊松过程

实际中事件强度常随时间变化，例如打车订单在夜间下降、白天上升。非齐次泊松过程（non-homogeneous Poisson process）用强度函数 $\lambda (t)\ge 0$ 取代常数 $\lambda$，要求增量独立，且

$$N_{t+s}-N_s\sim \mathrm{Poisson}({\int }_s^{s+t}\lambda (u)du)$$

事件在 $[0,T]$ 内发生的总数服从均值为 $\mathrm{\Lambda }(T)={\int }_0^T\lambda (u)du$ 的泊松分布。非齐次泊松过程保留了独立增量，放弃了平稳增量。若把时间按累积强度 $t\mapsto \mathrm{\Lambda }(t)$ 变换，过程化为强度为 1 的齐次泊松过程，这一变换称为时间变换，它把非齐次问题归约到齐次问题，也给出模拟非齐次泊松过程的一种算法（稀疏法）。

复合泊松过程与更新过程

复合泊松过程（compound Poisson process）在每次事件发生时附加一个随机标记。若事件到达构成强度为 $\lambda$ 的泊松过程，第 $k$ 次事件的标记为 $M_k$，且标记独立同分布、与到达时刻独立，则

$$S_t=\sum _{k=1}^{N_t}{M}_k$$

构成复合泊松过程。保险理赔总额、网络流量的累积字节数、某店铺单日订单总额都可以这样建模。由条件期望可以算出 $E[S_t]=\lambda tE[M_1]$，方差为 $\mathrm{Var}(S_t)=\lambda tE[M_1^2]$。

更新过程（renewal process）把到达间隔放宽为独立同分布的任意非负随机变量，不要求指数分布。更新函数 $m(t)=E[N_t]$ 满足更新方程

$$m(t)=F(t)+{\int }_0^{t}m(t-s)dF(s)$$

其中 $F$ 是间隔的分布函数。更新定理说明当 $t\to \mathrm{\infty }$ 时 $m(t)/t\to 1/\mu$，$\mu$ 为平均间隔，长期的平均事件率趋于平均间隔的倒数。泊松过程是指数间隔下的更新过程，更新过程把这一结构推广到一般间隔分布。

泊松过程与二项计数不能混用

在固定长度区间上统计事件数时，泊松模型与二项模型容易混淆。泊松模型描述的事件发生在连续时间上，事件数无上界，均值为 $\lambda t$。二项模型描述固定次数试验中的成功次数，试验次数事先固定。把单位时间内的事件数直接当作二项分布处理，会错误地引入一个并不存在的试验次数上限。判断依据是看事件能否在任意时刻发生：能则用泊松，试验次数固定则用二项。

import numpy as np

rng = np.random.default_rng(1)

# 齐次泊松过程：强度 lambda = 2，观察区间 [0, 10]
lam, T = 2.0, 10.0

# 方法一：按指数间隔累加到达时刻
arrivals = []
t = 0.0
while True:
    t += rng.exponential(1 / lam)
    if t >= T:
        break
    arrivals.append(t)
arrivals = np.array(arrivals)
print(f"区间内到达次数 N(10) = {len(arrivals)}")   # 理论期望 20

# 方法二：先抽总数 N(T) ~ Poisson(lambda * T)，再取均匀到达时刻
n = rng.poisson(lam * T)
times = np.sort(rng.uniform(0, T, size=n))
print(f"方法二的到达次数 = {n}")

# 验证到达间隔均值为 1/lambda
n_rep = 5000
mean_gaps = [rng.exponential(1 / lam, size=50).mean() for _ in range(n_rep)]
print(f"平均间隔的模拟均值 = {np.mean(mean_gaps):.4f}，理论值 = {1 / lam:.4f}")
# 0.5001 附近

# 验证 N(t) ~ Poisson(lambda t)：t = 3 时均值与方差均应为 6
counts = rng.poisson(lam * 3, size=100000)
print(f"N(3) 的模拟均值 = {counts.mean():.4f}，模拟方差 = {counts.var():.4f}，理论值均为 6")

2.8.5 布朗运动与相关过程

布朗运动的定义与性质

布朗运动（Brownian motion）$\{W_t{\}}_{t\ge 0}$ 是状态空间与指标集都连续的随机过程，由四条性质定义：$W_0=0$；增量独立，任意不重叠区间上的增量相互独立；增量正态，$W_t-W_s\sim \text{N}(0,t-s)$（$s<t$）；样本轨道连续。由增量正态可推出 $W_t\sim \text{N}(0,t)$，方差随时间线性增长。

布朗运动是随机游走在连续极限下的结果。把简单对称随机游走的时间与空间同时缩放，步长取 $1/\sqrt{n}$、步频取 $n$，轨道在 $n\to \mathrm{\infty }$ 时的极限就是布朗运动，这一结论称为唐斯克不变原理（Donsker's invariance principle）。它说明布朗运动的普适性，大量微小独立扰动叠加后的累积效应都趋向布朗运动。

布朗运动的轨道具有若干反直觉的性质。轨道处处连续，却在每一点都不可微，因此不能用普通微积分处理；轨道的二次变差在极限下等于时间长度，$[W{]}_t=t$，这一性质导致普通微积分的链式法则失效，需要引入伊藤积分。布朗运动具有自相似性，把时间缩放 $c$ 倍、空间缩放 $\sqrt{c}$ 倍后轨道分布不变。

几何布朗运动

几何布朗运动（geometric Brownian motion）由随机微分方程

$$d{S}_t=\mu S_{t}dt+\sigma S_{t}d{W}_t$$

定义，其中 $\mu$ 为漂移率，$\sigma$ 为波动率。它把布朗运动的比例增量作为模型，避免价格取负值。对 $\log S_t$ 用伊藤公式展开可得显式解

$$S_t=S_0\exp [(\mu -\frac{{\sigma }^2}{2})t+\sigma W_t]$$

因此 $S_t$ 服从对数正态分布，$\log S_t$ 服从正态分布。几何布朗运动是布莱克-斯科尔斯期权定价模型的基础资产模型，它的假设（对数收益正态、波动率恒定）与实际市场存在偏差，后续发展出随机波动率与跳跃扩散等模型。种群规模的变化也常用几何布朗运动描述，其中 $\mu$ 表示平均增长率，$\sigma$ 表示增长的随机波动。

分数布朗运动与长记忆性

分数布朗运动（fractional Brownian motion）$\{B_t^H\}$ 由 Hurst 指数 $H\in (0,1)$ 参数化，性质包括 $B_0^H=0$、增量平稳且高斯、自协方差函数为

$$\mathrm{Cov}(B_s^H,B_t^H)=\frac{1}{2}(s^{2H}+t^{2H}-|t-s{|}^{2H})$$

$H=1/2$ 时退化为布朗运动。$H>1/2$ 时增量的自相关为正且随间隔按幂律缓慢衰减，这一现象称为长记忆性；$H<1/2$ 时增量呈负相关，轨道比布朗运动更加曲折。$H\ne 1/2$ 时分数布朗运动的增量不独立，因此它不是半鞅，普通伊藤积分不适用。长记忆性在金融波动率、网络流量与生理信号（如心率变异性）中有实际应用。

2.8.6 高斯过程

定义

高斯过程（Gaussian process）是一族随机变量 $\{f(x){\}}_{x\in \mathcal{X}}$，其中指标集 $\mathcal{X}$ 可以是任意集合，要求任意有限个点 $x_1,\dots ,x_n$ 上的取值 $(f(x_1),\dots ,f(x_n))$ 服从联合正态分布。高斯过程由均值函数 $m(x)=E[f(x)]$ 与协方差函数 $k(x,x^{\prime })=\mathrm{Cov}(f(x),f(x^{\prime }))$ 完全确定，记作

$$f\sim \mathcal{GP}(m(\cdot ),k(\cdot ,\cdot ))$$

任意有限维边缘为正态这一要求把无穷维对象压缩为两个普通函数，这是高斯过程便于处理的根源。协方差函数 $k$ 必须满足对称性与半正定性：$k(x,x^{\prime })=k(x^{\prime },x)$，且任意点集上的 $n\times n$ 矩阵 $(k(x_i,x_j))$ 半正定。

核函数与协方差函数

协方差函数也称核函数（kernel function），它编码函数的先验性质，如光滑程度与周期结构。常用的核包括：

平方指数核（squared exponential kernel）$k(x,x^{\prime })={\sigma }^2\exp [-\frac{(x-x^{\prime }{)}^2}{2{\ell }^2}]$，其中 $\ell$ 为长度尺度，控制函数的起伏快慢，${\sigma }^2$ 为幅度。它对应的样本轨道无穷次可微，是最常用的默认核。

马特恩核（Matérn kernel）通过参数 $\nu$ 控制光滑度，$\nu =1/2$ 对应指数核，轨道粗糙；$\nu \to \mathrm{\infty }$ 时趋于平方指数核。空间统计与气象建模中常用 $\nu =3/2$ 或 $\nu =5/2$。

周期核（periodic kernel）刻画重复结构，适合建模昼夜节律或季节效应。把多个核相加或相乘还可以组合出趋势与噪声叠加的复合结构。

高斯过程回归

给定训练数据 $\{(x_i,y_i){\}}_{i=1}^n$，设观测为函数值加噪声 $y_i=f(x_i)+{\epsilon }_i$，${\epsilon }_i\sim \text{N}(0,{\sigma }_n^2)$。由高斯过程的定义，训练值与测试点 $x_{\ast }$ 上的函数值联合服从正态分布：

$$\left[\begin{array}{c}\mathbf{y}\\ f_{\ast }\end{array}\right]\sim \text{N}(\mathbf{0},\left[\begin{array}{cc}K+{\sigma }_n^{2}I& {\mathbf{k}}_{\ast }\\ {\mathbf{k}}_{\ast }^{\mathrm{\top }}& k_{\ast \ast }\end{array}\right])$$

其中 $K$ 是训练点上的核矩阵，${\mathbf{k}}_{\ast }=(k(x_{\ast },x_i))$，$k_{\ast \ast }=k(x_{\ast },x_{\ast })$。由正态分布的条件分布公式得到后验预测分布

$$f_{\ast }\mid \mathbf{y}\sim \text{N}({\mathbf{k}}_{\ast }^{\mathrm{\top }}(K+{\sigma }_n^{2}I{)}^{-1}\mathbf{y},k_{\ast \ast }-{\mathbf{k}}_{\ast }^{\mathrm{\top }}(K+{\sigma }_n^{2}I{)}^{-1}{\mathbf{k}}_{\ast })$$

均值是训练值的线性组合，方差随测试点远离训练点而增大，在无数据处回到先验方差。这一公式把回归问题写成了闭式解，代价是求 $n\times n$ 矩阵的逆，复杂度为 $O(n^3)$。超参数 $\ell ,{\sigma }^2,{\sigma }_n^2$ 通过最大化对数边缘似然确定，边缘似然自动包含对模型复杂度的惩罚，无需额外的正则化项。

高斯过程分类

分类任务中观测是离散标签，高斯过程不能直接给出标签分布。做法是对每个类别引入一个潜在函数 $f$，通过链接函数把 $f$ 映射为类别概率，二分类时取 $\pi (x)=\mathrm{\Phi }(f(x))$ 或 Sigmoid 函数。后验分布不再有闭式解，需要用拉普拉斯近似或期望传播等近似推断方法。预测时对潜在函数求积分得到类别概率。高斯过程分类在样本量较小的场景中表现稳定，计算复杂度同样为 $O(n^3)$，样本量很大时需要用稀疏近似。

与贝叶斯非参数的联系

高斯过程是贝叶斯非参数（Bayesian nonparametrics）的典型例子。参数模型的参数个数固定，高斯过程把先验放在函数空间上，模型的有效复杂度随数据量自动增长。它与再生核希尔伯特空间（reproducing kernel Hilbert space）有紧密联系：给定核 $k$，存在唯一的希尔伯特空间，高斯过程回归的后验均值落在该空间中。核岭回归、支持向量机与高斯过程回归都可以在这一框架下理解，它们使用相同的核，区别在于是否给出完整的后验分布。

import numpy as np

rng = np.random.default_rng(3)

def rbf(a, b, length=1.0, sigma=1.0):
    """平方指数核"""
    d = a[:, None] - b[None, :]
    return sigma**2 * np.exp(-0.5 * (d / length) ** 2)

# 训练数据
X = np.array([-3.0, -1.5, 0.0, 1.5, 3.0])
y = np.sin(X) + 0.1 * rng.standard_normal(len(X))
noise = 0.01

K = rbf(X, X) + noise * np.eye(len(X))
Ks = np.linalg.solve(K, y)

# 测试点上的后验均值与方差
Xstar = np.linspace(-4, 4, 9)
Kss = rbf(Xstar, X)
mu = Kss @ Ks
cov = rbf(Xstar, Xstar) - Kss @ np.linalg.solve(K, Kss.T)
print(f"后验均值 = {np.round(mu, 3)}")
print(f"后验标准差 = {np.round(np.sqrt(np.diag(cov)), 3)}")
# 靠近训练点处标准差小，远离处趋向先验标准差 1.0

2.8.7 鞅与可选停止定理

鞅的定义与例子

鞅（martingale）是一类增量均值为零的过程。设 $\{X_n\}$ 是随机过程，$\{{\mathcal{F}}_n\}$ 是递增的信息流（$\sigma$-代数），且 $X_n$ 关于 ${\mathcal{F}}_n$ 可测、可积，若对任意 $n$，

$$E[X_{n+1}\mid {\mathcal{F}}_n]=X_n$$

则称 $\{X_n\}$ 为鞅。把等号改为 $\ge$ 得到下鞅，改为 $\le$ 得到上鞅。信息流 ${\mathcal{F}}_n$ 表示到时刻 $n$ 为止可获得的全部信息，条件期望在第 $n$ 步信息下给出对下一步的最优预测。鞅的含义是当前值是下一步的期望，过程没有系统性上升或下降的趋势。

简单对称随机游走是鞅：$E[S_{n+1}\mid S_n]=S_n$，因为下一步的增量为 $+1$ 与 $-1$ 的概率各半。公平赌博中赌徒的财富是鞅，因为每次下注的期望收益为零。若赌局对赌徒有利则财富是下鞅，不利则为上鞅。折扣后的资产价格在风险中性测度下是鞅，这是金融衍生品定价的核心。平方可积的鞅序列 $\{S_n^2-n\}$ 也是鞅，因为它减去累积方差后期望不再增长，这一构造是计算随机游走首达时分布的工具。

可选停止定理

可选停止定理（optional stopping theorem）把鞅的等期望性质从固定时刻推广到随机时刻。设 $\{X_n\}$ 是鞅，$\tau$ 是停时（即 $\{\tau =n\}\in {\mathcal{F}}_n$，是否停止在时刻 $n$ 可由当前信息判定），若满足下列条件之一，则 $E[X_{\tau }]=E[X_0]$：$\tau$ 有界；$\tau$ 几乎必然有限且增量有界；$\{X_n\}$ 一致可积。定理的意义是把公平游戏在任意合理停止策略下平均收益仍为零这一直觉严格化。

定理的条件不能省略。赌徒停止策略的例子说明违背条件时结论失效：在简单对称随机游走中，采用首次达到 $+1$ 就停止的策略，$\tau$ 有限，且 $E[S_{\tau }]=1\ne 0=E[S_0]$，原因是 $\tau$ 无界且停止时的取值无界，可选停止定理的条件不成立。这一反例常被误读为存在必胜策略，实际原因在于策略需要无限资金与无限时间，一旦加上资金与时间上限，等期望性质恢复。

2.8.8 平稳过程与谱分析

严平稳与宽平稳

严平稳（strictly stationary）要求过程的有限维分布在时间平移下不变：对任意 $n$，$t_1,\dots ,t_n$ 与任意 $h$，$(X_{t_1},\dots ,X_{t_n})$ 与 $(X_{t_1+h},\dots ,X_{t_n+h})$ 同分布。严平稳条件很强，实际验证困难。

宽平稳（weakly stationary）只要求前两阶矩满足条件：均值函数为常数 $E[X_t]=\mu$；自协方差只依赖时间差，$C(s,t)=\gamma (t-s)$。对高斯过程而言，宽平稳与严平稳等价，因为正态分布完全由前两阶矩决定。宽平稳的条件容易验证，也足够支撑自相关与谱分析，实践中的平稳性通常指宽平稳。若过程经差分后平稳，称其为差分平稳，这是单位根检验的基础。

各态历经性

各态历经性（ergodicity）指时间平均等于空间平均。对宽平稳过程，若均值各态历经成立，则单条长轨道的时间平均依概率收敛到总体均值：

$$\frac{1}{T}{\int }_0^T{X}_{t}dt⟶\mu$$

各态历经性使我们可以从一次观测推知总体矩，这是时间序列分析的立足点。均值各态历经的一个充分条件是自协方差函数可积，${\int }_0^{\mathrm{\infty }}|\gamma (h)|dh<\mathrm{\infty }$，它要求远处的相关性足够快地衰减。

自相关与偏自相关函数

宽平稳过程的自协方差函数 $\gamma (h)=\mathrm{Cov}(X_t,X_{t+h})$，自相关函数（autocorrelation function，ACF）定义为

$$\rho (h)=\frac{\gamma (h)}{\gamma (0)}$$

它给出相隔 $h$ 步的取值之间的相关系数，$|\rho (h)|\le 1$，$\rho (0)=1$。ACF 刻画过程的记忆长度与周期结构，衰减快的 ACF 对应短记忆。

偏自相关函数（partial autocorrelation function，PACF）度量在去掉中间变量影响后 $X_t$ 与 $X_{t+h}$ 的相关性。$h$ 阶偏自相关是回归 $X_t$ 对 $X_{t-1},\dots ,X_{t-h+1}$ 之后 $X_{t-h}$ 的系数。AR 模型的 ACF 缓慢衰减、PACF 在阶数处截断，MA 模型的 ACF 截断、PACF 缓慢衰减，这一对照是识别模型阶数的常用手段。

功率谱密度

功率谱密度（power spectral density）把过程按频率分解。由维纳-辛钦定理（Wiener-Khinchin theorem），宽平稳过程的自协方差与谱密度构成傅里叶变换对：

$$S(\omega )=\sum _{h=-\mathrm{\infty }}^{\mathrm{\infty }}\gamma (h)e^{-i\omega h},\gamma (h)=\frac{1}{2\pi }{\int }_{-\pi }^{\pi }S(\omega )e^{i\omega h}d\omega$$

谱密度 $S(\omega )$ 描述各频率成分的功率大小。周期性成分在谱上表现为尖峰，据此可以识别心跳、呼吸或机器振动中的周期信号。谱分析把时域的卷积关系化为频域的乘积关系，线性滤波器的分析因此变得简单。

白噪声与有色噪声

白噪声（white noise）是均值为零、方差恒定且不同时刻不相关的序列，$\gamma (h)=0$（$h\ne 0$）。由维纳-辛钦定理，白噪声的谱密度为常数，各频率的功率相等。高斯白噪声是最常见的噪声模型，也是时间序列建模中残差的标准假设。有色噪声（colored noise）指不同时刻相关的噪声，其谱密度随频率变化。1/f 噪声在低频处功率较大，出现在心率变异性与神经信号中；有色噪声常被建模为 MA 或 ARMA 过程，通过预白化把有色噪声转化为白噪声后再做检测。

2.8.9 概率图模型

贝叶斯网络

概率图模型（probabilistic graphical model）用图结构编码随机变量之间的条件独立关系。贝叶斯网络（Bayesian network）是有向无环图（directed acyclic graph，DAG），节点对应随机变量，有向边表示直接依赖。给定父节点后，每个节点与其非后代节点条件独立，这一性质使联合分布分解为

$$P(X_1,\dots ,X_n)=\prod _{i=1}^{n}P(X_i\mid \mathrm{Pa}(X_i))$$

其中 $\mathrm{Pa}(X_i)$ 是节点 $i$ 的父节点集合。分解把联合分布的参数量从指数级降为各条件分布参数量之和。以链式法则为起点，贝叶斯网络通过条件独立的假设删去了多余的因子，这与 2.1 节的链式法则一脉相承。节点之间的条件独立关系可以用 d-分离判定：若一条路径上存在被观测节点阻断的对撞结构（两个箭头指向同一节点且该节点及其后代未被观测），则路径被阻断。

马尔可夫随机场与条件随机场

马尔可夫随机场（Markov random field，MRF）用无向图表示依赖关系，适合描述对称的相互作用，例如图像中相邻像素的相似性。无向图不能像有向图那样为每个节点指定父节点，因子分解改用团（clique，两两相连的节点集合）上的势函数（potential function） ${\psi }_C$：

$$P(\mathbf{x})=\frac{1}{Z}\prod _{C\in \mathcal{C}}{\psi }_C({\mathbf{x}}_C),Z=\sum _{\mathbf{x}}\prod _{C\in \mathcal{C}}{\psi }_C({\mathbf{x}}_C)$$

归一化常数 $Z$ 称为配分函数，需对所有变量组合求和，计算量随变量数指数增长，这一困难是无向图模型推断的核心难点。势函数只要求非负，不需要像概率那样归一化，全局归一化由 $Z$ 一次性完成。

条件随机场（conditional random field，CRF）是判别式模型，直接建模条件分布

$$P(\mathbf{y}\mid \mathbf{x})=\frac{1}{Z(\mathbf{x})}\prod _C{\psi }_C({\mathbf{y}}_C,\mathbf{x})$$

势函数依赖观测 $\mathbf{x}$，归一化只对标签 $\mathbf{y}$ 进行。CRF 不需要对观测分布建模，放宽了 HMM 中观测条件独立的假设，序列标注任务中常取得更好的效果。线性链 CRF 在给定观测后结构退化为一条链，可用前向-后向算法高效计算，这与 HMM 的算法结构相同。

因子图与信念传播

因子图（factor graph）是二分图，同时包含变量节点与因子节点，因子节点连到它涉及的变量。因子图统一了有向图与无向图的表示，任何概率图模型都可以转化为因子图。在因子图上，信念传播（belief propagation）算法通过节点之间传递消息来近似计算边缘分布。

消息从变量节点传向因子节点或反向，每条消息是一个关于某个变量的函数。变量节点对所有入消息取乘积并归一化，因子节点对涉及的其他变量求和（或积分）后传给目标变量。对于树结构的因子图，信念传播给出精确的边缘分布；对于含环的图，迭代传播给出近似结果，称为循环信念传播。和-积算法计算边缘分布，最大-积算法计算最大后验配置，维特比算法是最大-积算法在链结构上的特例。

变分消息传递

图结构复杂时，精确推断的计算量无法承受。变分消息传递（variational message passing）用变分推断近似后验。基本思路是在一个受限制的分布族（例如各变量相互独立的平均场族）中寻找与真实后验最接近的分布，以最小化库尔贝克-莱布勒散度为准则。更新规则表现为节点之间传递消息，形式与信念传播相似，含义是把信念传播中的求和替换为对近似后验的期望。变分方法把推断问题转化为优化问题，通常比精确推断快，代价是引入近似误差。变分推断在 2.12 节与变分自编码器中都会用到。

图模型是共同语言

概率图模型的价值在于把不同模型的共同结构抽象出来。HMM 是链结构的贝叶斯网络，卡尔曼滤波是线性高斯状态空间模型，条件随机场是判别式的无向链，主题模型是含隐变量的有向图。一旦用图的语言表达，推断算法、参数估计算法与近似方法就可以在这些模型之间迁移。掌握因子图与消息传递，等于掌握了一类推断问题的通用解法。

2.8.10 因果图与结构因果模型简介

有向无环图与因果结构

结构因果模型（structural causal model，SCM）用一组结构方程描述变量之间的因果关系。每个变量由它的直接原因与一个独立的外生噪声决定：

$$X_i=f_i(\mathrm{Pa}(X_i),U_i)$$

其中 $U_i$ 是相互独立的噪声。结构方程的方向具有因果含义，箭头从原因指向结果，与贝叶斯网络中箭头只表示条件依赖不同。变量之间的因果结构用有向无环图表示，图中不存在路径从祖先绕回自身。图结构决定了哪些变量可以作为控制变量：控制对撞节点会引入虚假关联，控制中介节点会阻断部分因果效应，这些规则由 d-分离与后门准则给出。

干预与 do 算子的直观含义

观测与干预是两类不同的操作。观测 $X=x$ 只是筛选出 $X$ 取该值的样本，其余变量仍按原分布变化；干预 $do(X=x)$ 把 $X$ 的取值固定，切断 $X$ 与它原有原因之间的联系，让其余变量的分布按结构方程重新生成。因果效应用干预分布 $P(Y\mid do(X=x))$ 刻画，它回答的是改变 $X$ 之后 $Y$ 会怎样这一问题，与观测条件分布 $P(Y\mid X=x)$ 一般不同。

do 算子把干预形式化。在因果图中，$do(X=x)$ 相当于删除所有指向 $X$ 的边，把 $X$ 固定为 $x$，再在修改后的图上计算 $Y$ 的分布。当存在一组变量满足后门准则（阻断了所有从 $X$ 到 $Y$ 的后门路径）时，干预分布可以通过对这些变量调整后从观测数据估计：

$$P(Y\mid do(X=x))=\sum _{z}P(Y\mid X=x,Z=z)P(z)$$

公式称为调整公式，它是观察数据推断因果效应的桥梁。随机对照试验通过随机化切断混杂变量的影响，相当于直接执行 $do$ 操作，因此是估计因果效应的金标准。

因果推断的完整内容，包括后门准则、前门准则、工具变量与反事实推理，将在 2.14 节展开。本节的作用是说明因果图与概率图模型共享同一套图语言：节点、边、条件独立、因子分解与消息传递都可以迁移到因果结构的分析中。贝叶斯网络刻画变量之间的依赖与独立关系，因果图在此基础上增加方向与干预的语义，两者共同构成后续因果推断与贝叶斯建模的语言基础。

2.8.11 本节小结

从概率分布到结构化的随机过程

本节把随机变量的视角从有限维扩展到指标集上的族。随机过程由有限维分布族确定，均值与自协方差刻画其二阶结构。马尔可夫链用无后效性把依赖压缩为转移矩阵，$n$ 步转移满足查普曼-科尔莫戈罗夫方程，状态的常返、周期与互通类决定长期行为，平稳分布与细致平衡给出不变性，遍历性保证时间平均趋于空间平均，这两条性质分别是 PageRank 与 MCMC 的理论基础。隐马尔可夫模型把状态链与观测层分离，前向算法计算似然，后向算法与维特比算法完成平滑与解码，鲍姆-韦尔奇算法作为 EM 的特例完成参数估计。泊松过程用独立平稳增量刻画随机到达，到达间隔服从指数分布，非齐次与复合形式拓宽了适用场景，更新过程把间隔分布放宽为一般分布。布朗运动是随机游走的连续极限，几何布朗运动给出资产价格模型，分数布朗运动刻画长记忆。高斯过程把先验放在函数空间上，核函数决定函数性质，后验预测有闭式解，它是贝叶斯非参数的入口。鞅与可选停止定理描述无偏随机过程的停止性质，平稳过程与谱分析把时域相关性转化为频域功率，白噪声是其中的基准。概率图模型用有向图与无向图统一表达条件独立与因子分解，因子图上的信念传播与变分消息传递提供推断算法。因果图在依赖结构之上增加干预语义，为后续因果推断提供共同语言。

下一节转向统计推断。前面的内容都在给定分布或过程的条件下讨论性质，而实际建模时分布本身往往未知，需要从观测数据推断。2.9 节给出统计推断的基本设定，讨论抽样分布、估计量的评价标准、充分性、Fisher 信息与 Cramér-Rao 下界，以及刀切法与自助法这类重抽样工具，为点估计与区间估计建立框架。

练习题

第 1 题 概念推导

设 $\{X_n{\}}_{n\ge 0}$ 是状态空间为 $S$ 的时齐马尔可夫链，转移矩阵为 $P$。证明 $n$ 步转移概率满足查普曼-科尔莫戈罗夫方程 $p_{ij}^{(m+n)}=\sum _{k\in S}{p}_{ik}^{(m)}{p}_{kj}^{(n)}$，并由此说明为什么 $n$ 步转移矩阵等于 $P^n$。

参考答案

用全概率公式按时刻 $m$ 的状态展开。在时刻 $0$ 处于状态 $i$ 的条件下，时刻 $m$ 必然处于某个状态 $k\in S$，且这些事件构成一个划分。由马尔可夫性质，从 $i$ 到 $k$ 的 $m$ 步转移只依赖 $i$ 与 $k$，从 $k$ 到 $j$ 的 $n$ 步转移只依赖 $k$ 与 $j$，二者条件独立。因此

$$p_{ij}^{(m+n)}=P(X_{m+n}=j\mid X_0=i)=\sum _{k\in S}P(X_m=k\mid X_0=i)P(X_{m+n}=j\mid X_m=k)=\sum _{k\in S}{p}_{ik}^{(m)}{p}_{kj}^{(n)}$$

右边正是矩阵乘积 $P^m{P}^n$ 的第 $(i,j)$ 项。由 $m,n$ 的任意性，方程对一切非负整数成立，取 $m=1$、$n=1$ 得 $P^{(2)}=P^2$，归纳即得 $P^{(n)}=P^n$。方程把多步转移的计算归结为矩阵乘法，链的长期行为因此可以从单个一步转移矩阵读取。

第 2 题 计算推理

某马尔可夫链有两个状态 $0,1$，转移矩阵为

$$P=\left[\begin{array}{cc}1-\alpha & \alpha \\ \beta & 1-\beta \end{array}\right],0<\alpha ,\beta <1$$

求平稳分布 $\pi =({\pi }_0,{\pi }_1)$，判断链是否可逆，并说明当 $\alpha =\beta$ 时平稳分布的取值。

参考答案

平稳分布满足 $\pi P=\pi$ 与 ${\pi }_0+{\pi }_1=1$。展开第一行：

$${\pi }_0(1-\alpha )+{\pi }_1\beta ={\pi }_0$$

整理得 ${\pi }_1\beta ={\pi }_0\alpha$，即 ${\pi }_0\alpha ={\pi }_1\beta$。结合归一化条件 ${\pi }_0+{\pi }_1=1$ 解得

$${\pi }_0=\frac{\beta }{\alpha +\beta },{\pi }_1=\frac{\alpha }{\alpha +\beta }$$

验证细致平衡：${\pi }_0{P}_{01}={\displaystyle \frac{\beta }{\alpha +\beta }}\alpha ={\displaystyle \frac{\alpha \beta }{\alpha +\beta }}$，${\pi }_1{P}_{10}={\displaystyle \frac{\alpha }{\alpha +\beta }}\beta ={\displaystyle \frac{\alpha \beta }{\alpha +\beta }}$，两者相等，因此链是可逆的。当 $\alpha =\beta$ 时，${\pi }_0={\pi }_1=1/2$，平稳分布为均匀分布，链在长期停留于两个状态的时间各占一半，此时转移矩阵是对称的。链不可约且非周期，平稳分布即极限分布，与初始状态无关。

第 3 题 代码验证

用蒙特卡洛方法验证布朗运动的二次变差性质。取长度为 $1$ 的时间区间，把它分为 $n$ 个小区间，生成布朗运动的离散路径，计算增量平方和 $\sum _i(W_{t_{i+1}}-W_{t_i}{)}^2$。分别在 $n=100,1000,10000$ 下重复多次，观察增量平方和的均值如何随 $n$ 变化，并与理论值 $[W{]}_1=1$ 对照。

参考答案

import numpy as np

rng = np.random.default_rng(42)
T = 1.0
n_rep = 2000

for n in [100, 1000, 10000]:
    quad = []
    for _ in range(n_rep):
        # 布朗运动在等分点上的取值：增量独立，方差为 dt
        dt = T / n
        inc = rng.normal(0, np.sqrt(dt), size=n)
        quad.append(np.sum(inc ** 2))
    print(f"n = {n:6d}，增量平方和均值 = {np.mean(quad):.4f}，标准差 = {np.std(quad):.4f}")
# n =    100，增量平方和均值 = 1.0003，标准差 = 0.1414
# n =   1000，增量平方和均值 = 1.0005，标准差 = 0.0448
# n =  10000，增量平方和均值 = 1.0002，标准差 = 0.0141

三种划分下增量平方和的均值都接近 $1$，与理论值 $[W{]}_1=1$ 一致，说明布朗运动是有限二次变差过程。标准差随 $n$ 增大按 $1/\sqrt{n}$ 衰减，因为平方和是 $n$ 个独立同分布项的累加。这一性质说明布朗运动的轨道虽然连续，但其增量平方不会随时间消失，普通微积分的可微性假设不成立，随机积分需要单独建立。

第 4 题 综合应用

某网约车平台把订单到达建模为泊松过程，平均每小时到达 $6$ 单。设工作日下午时段为 $14:00$ 到 $18:00$，强度在此时段恒为 $6$，夜间 $22:00$ 到次日 $06:00$ 强度降为 $2$。计算下列量：（一）下午时段内到达订单数的均值与标准差；（二）下午时段内两次到达间隔的均值；（三）夜间时段到达订单数的均值。若观察到某天下午到达订单数为 $40$ 单，说明如何据此检验泊松假设。

参考答案

（一） 下午时段长度为 $4$ 小时，强度 $\lambda =6$ 单/小时，到达订单数 $N$ 服从泊松分布，参数 $\mathrm{\Lambda }=\lambda t=24$，因此

$$E[N]=24,\mathrm{Std}(N)=\sqrt{24}\approx 4.90$$

（二） 到达间隔服从参数为 $\lambda =6$ 的指数分布，均值为 $1/\lambda =1/6$ 小时，即约 $10$ 分钟。

（三） 夜间时段长度为 $8$ 小时，强度降为 $\lambda =2$，参数 $\mathrm{\Lambda }=2\times 8=16$，到达订单数均值为 $16$。

检验泊松假设：泊松分布的一个重要特征是均值与方差相等。若某天下午到达 $40$ 单，在均值为 $24$ 的泊松假设下，$40$ 偏离均值约 $(40-24)/\sqrt{24}\approx 3.27$ 个标准差，属于较罕见的事件。收集多天的到达数量后，可以比较样本均值与样本方差，若方差明显大于均值（过离散），说明到达并非完全随机，可能受到未建模的时段效应影响，此时应改用非齐次泊松过程或引入过离散模型。单个观测不足以判断，需要累积足够样本再作结论。

常见错误

错误 1 · 把平稳分布当作任意链都存在且唯一

原因：只有不可约、正常返的链才存在唯一的平稳分布。可约链的每个闭合互通类上各有一个平稳分布，存在无穷多个平稳分布的凸组合；零常返链没有平稳概率分布。把方程 $\pi =\pi P$ 有解直接等同于平稳分布存在且唯一，会在可约链上给出错误结论。

解决：先判断链的互通类结构与常返性质，再讨论平稳分布。有限状态的链只需要检查不可约性即可保证存在唯一平稳分布，因为有限不可约链必为正常返。解平稳分布时若方程有多组线性无关解，说明链可约，需要补充实际关心的初始条件下链落入哪个类。

错误 2 · 混淆观测条件分布与干预分布

原因：$P(Y\mid X=x)$ 与 $P(Y\mid do(X=x))$ 是两个不同的量，前者包含混杂变量造成的虚假关联，后者反映改变 $X$ 的因果效应。在因果图中存在后门路径时，两者不相等。把回归系数直接解释为因果效应，属于这一混淆的常见形式。

解决：先画出因果图，找出从 $X$ 到 $Y$ 的路径，确定哪些是前门路径、哪些是后门路径。若存在后门路径，需要选取满足后门准则的变量集 $Z$ 并调整，用 $\sum _{z}P(Y\mid X=x,Z=z)P(z)$ 估计干预分布。没有因果图支持时，回归系数只能作关联解释。

错误 3 · 在非平稳序列上直接套用平稳过程工具

原因：自相关函数、功率谱密度与遍历定理都以宽平稳为前提。存在趋势或单位根的序列不是平稳序列，直接计算 ACF 会得到缓慢衰减的假象，谱分析会给出低频处的虚假功率，时间平均也不再收敛到总体均值。

解决：分析前先检查平稳性，画出序列图观察趋势与方差变化，用单位根检验判断差分次数。对确定性趋势可以做去趋势处理，对随机趋势做差分，直到序列平稳后再计算 ACF、PACF 与谱密度。建模时优先使用能容纳非平稳性的模型。

错误 4 · 忽略 HMM 的观测条件独立假设

原因：HMM 假设给定隐状态后各时刻观测相互独立，发射概率只依赖当前状态。当观测存在重叠特征或长距离相关时，这一假设明显违背，前向算法给出的似然与维特比解码结果会产生偏差，模型输出的概率往往过于自信。

解决：需要保留马尔可夫结构又想放宽观测假设时，改用条件随机场直接建模 $P(Y\mid X)$，它通过全局归一化避免对观测分布建模。也可以扩展隐状态空间，把观测中相关的部分纳入状态，或用重叠特征构造发射分布。选择模型时先评估观测独立假设是否合理，再决定是否采用 HMM。