随机过程深度实战：从随机游走、泊松过程与布朗运动到伊藤积分、马尔可夫链与高斯过程的完整工程链路

随机过程是把"不确定性"变成可计算对象的第一性原理工具。我们日常写的神经网络、队列系统、扩散模型、强化学习环境、贝叶斯先验，背后几乎都站着某个随机过程——SGD 的噪声是布朗扰动，MCMC 采样是马尔可夫链，扩散模型的反向去噪是伊藤随机微分方程，PageRank 是平稳分布，主动学习的采集函数是高斯过程。本文用纯 NumPy 从零实现随机游走、泊松过程、布朗运动、伊藤 SDE、马尔可夫链、鞅与高斯过程七条主线，把抽象定义逐个翻译成可运行、可验证的工程代码，并在结尾给出与 AI 系统的完整映射表。

一、第一性原理：什么是随机过程

一个随机过程是一族随机变量 {X(t), t ∈ T}，其中 T 是指标集（时间，离散或连续），每个 X(t) 取值于状态空间 S（离散或连续）。一次具体观测得到的是一条样本路径 ω → {X(t, ω)}。

工程上最关心的三类结构：

• 增量结构：独立增量（泊松、布朗）、马尔可夫增量（只依赖当前状态）、鞅增量（条件期望为 0）。
• 滤波（filtration） F_t：到时刻 t 为止所有可见信息的累积，是"条件期望"的严格载体。
• 尺度变换：布朗运动的 W(ct) ≡ c^{1/2} W(t)、泊松过程的叠加/稀释，决定了过程如何在大时间尺度上重标定。

判定一个过程"可工程化"的关键是：它有没有封闭的分布刻画（随机游走的 CLT、泊松的计数律、布朗的二次变差、马尔可夫的平稳分布、GP 的高斯后验）。下面逐一落地。

二、随机游走与中心极限定理

最简单的随机过程是对称随机游走 S_n = Σ ξ_i，ξ_i ∈ {±1} 等概率独立。它的均值恒为 0、方差 n，而标准化后 S_n/√n 收敛到标准正态分布——这正是布朗运动的离散起源。

import numpy as np

# 1) 中心极限定理：S_n / sqrt(n) -> N(0,1)
N = 2000
M = 50000
rng = np.random.default_rng(42)
S = rng.choice([-1, 1], size=(M, N)).cumsum(axis=1)
scaled = S[:, -1] / np.sqrt(N)
print("mean =", round(scaled.mean(), 4))            # 理论 0
print("std  =", round(scaled.std(ddof=1), 4))       # 理论 1
q = np.percentile(scaled, [2.5, 50, 97.5])
print("empirical 2.5/50/97.5% =", np.round(q, 3))  # 应接近 -1.96 / 0 / 1.96

# 2) 首达时：对称游走首次触及 ±a 的期望为 a^2（仅到达单侧 +a 的期望发散）
def first_passage(a, max_n=500000, seed=7):
    rng = np.random.default_rng(seed)
    s = 0
    for n in range(1, max_n + 1):
        s += rng.choice([-1, 1])
        if abs(s) >= a:          # 触及 +a 或 -a 即停
            return n
    return None

fpt = np.array([first_passage(20, seed=i) for i in range(3000)])
print("E[T_20] ~", round(fpt.mean(), 1), "(理论 a^2 = 400)")

输出应显示标准化序列均值≈0、标准差≈1、分位数逼近 ±1.96，且首达时均值逼近 400。CLT 告诉我们：即使每一步毫无方向，大样本下路径的整体形态仍然可预测——这是蒙特卡洛方法的根基。

三、泊松过程：计数与到达

泊松过程是"独立平稳增量 + 单位时间平均 λ 次事件"的连续时间计数过程。N(t) ~ Poisson(λt)，到达间隔 Δt ~ Exp(λ)，且两个独立泊松过程叠加仍是泊松过程（速率相加）。它几乎是"事件流"建模的默认选择：请求到达、异常爆发、光子计数。

import numpy as np

def poisson_process(lam, T, seed=0):
    rng = np.random.default_rng(seed)
    times, t = [], 0.0
    while t < T:
        t += rng.exponential(1.0 / lam)
        if t < T:
            times.append(t)
    return np.array(times)

lam, T = 3.0, 1000.0
counts = np.array([len(poisson_process(lam, T, seed=i)) for i in range(5000)])
print("mean N(T) =", round(counts.mean(), 2), "(理论 λT =", lam * T, ")")
print("var  N(T) =", round(counts.var(ddof=1), 2), "(理论 λT =", lam * T, ")")

arr = poisson_process(lam, T, seed=99)
inter = np.diff(arr)
print("mean interarrival =", round(inter.mean(), 4), "(理论 1/λ =", 1 / lam, ")")

# 叠加性：PP(1.5) + PP(1.5) -> PP(3.0)
a = poisson_process(1.5, T, seed=11)
b = poisson_process(1.5, T, seed=12)
merged = np.sort(np.concatenate([a, b]))
print("merged rate ~", round(len(merged) / T, 3), "(理论 3.0)")

N(T) 的均值与方差都应逼近 3000，到达间隔均值逼近 0.333，叠加后速率逼近 3.0。均值=方差正是泊松分布的特征指纹，也是区分"泊松噪声"与"过离散噪声"的工程判据。

四、布朗运动与二次变差（伊藤微积分的基石）

布朗运动 W(t) 满足：W(0)=0、路径连续、增量 W(t)-W(s) ~ N(0, t-s) 且独立平稳。它最反直觉的性质是二次变差 [W]_t = lim Σ (ΔW_i)² = t——平方增量不消失，反而累积成时间本身。这一条直接催生了伊藤微积分与普通微积分的根本分野。

import numpy as np

def brownian(n, dt=1e-4, seed=0):
    rng = np.random.default_rng(seed)
    dW = rng.normal(0, np.sqrt(dt), size=n)
    return np.cumsum(dW), dW

n, dt = 100000, 1e-4
T = n * dt
W, dW = brownian(n, dt, seed=3)
print("sum dW^2 =", round(np.sum(dW ** 2), 3), " 目标 T =", T)   # 二次变差 -> T

rng = np.random.default_rng(5)
WT = np.array([brownian(n, dt, seed=i)[0][-1] for i in range(2000)])
print("E[W_T^2] =", round((WT ** 2).mean(), 3), " 目标 T =", T)
print("E[W_T]   =", round(WT.mean(), 4))                         # 应≈0

Σ dW² 应逼近 10（即 T），E[W_T²] 同样逼近 10，而 E[W_T]≈0。二次变差是连接"噪声"与"时间"的桥梁：它让伊藤引理多出一个 dt 项，也让随机微分方程与常微分方程彻底分道扬镳。

五、伊藤积分与几何布朗运动（随机微分方程的落地）

伊藤微分把布朗增量当作"噪声源"：对 dX = μ dt + σ dW，普通链式求导失效，必须使用伊藤引理 df = (f_t + μ f_x + ½σ² f_xx) dt + σ f_x dW。几何布朗运动 dS = μ S dt + σ S dW 是金融与扩散模型的核心，其解为 S_t = S_0 exp((μ-σ²/2)t + σ W_t)——注意漂移被 -σ²/2 修正，这是伊藤修正项的直接体现。

import numpy as np

mu, sigma, S0, T, n = 0.1, 0.3, 100.0, 1.0, 5000

def gbm_euler(seed):
    rng = np.random.default_rng(seed)
    dt = T / n
    dW = rng.normal(0, np.sqrt(dt), size=n)
    S = S0
    for dW_i in dW:
        S = S + mu * S * dt + sigma * S * dW_i          # 伊藤欧拉-丸山
    return S

paths = np.array([gbm_euler(i) for i in range(4000)])
print("E[S_T]      =", round(paths.mean(), 2),
      " 理论 S0*exp(μT) =", round(S0 * np.exp(mu * T), 2))
print("E[log S_T]  =", round(np.log(paths).mean(), 4),
      " 理论 log S0 + (μ-σ²/2)T =", round(np.log(S0) + (mu - sigma ** 2 / 2) * T, 4))

E[S_T] 应逼近 110.5，E[log S_T] 应逼近 4.660。两者差了一个 σ²/2——若误用普通微分的"无修正"漂移，对数收益会系统性高估。这正是扩散模型（DDPM 的反向 SDE）与随机优化里必须区分伊藤/斯特拉托诺维奇的深层原因。

六、马尔可夫链：状态转移与平稳分布

马尔可夫性要求 P(X_{n+1} | X_n, …, X_0) = P(X_{n+1} | X_n)——未来只依赖现在。长期行为由平稳分布 π 决定，满足 π P = π。它同时是 PageRank、MCMC 采样与强化学习环境动力学的统一语言。

import numpy as np

# 反射边界上的生灭链（对称转移）
P = np.array([
    [0.5,  0.5,  0.0,  0.0,  0.0],
    [0.25, 0.5,  0.25, 0.0,  0.0],
    [0.0,  0.25, 0.5,  0.25, 0.0],
    [0.0,  0.0,  0.25, 0.5,  0.25],
    [0.0,  0.0,  0.0,  0.5,  0.5],
])

def stationary_power(P, iters=5000):
    pi = np.ones(P.shape[0]) / P.shape[0]
    for _ in range(iters):
        pi = pi @ P
    return pi

def stationary_lin(P):
    n = P.shape[0]
    A = P.T - np.eye(n)
    A[-1] = 1.0
    b = np.zeros(n); b[-1] = 1.0
    return np.linalg.solve(A, b)

pi1 = stationary_power(P)
pi2 = stationary_lin(P)
print("幂迭代 :", np.round(pi1, 4))
print("线性方程:", np.round(pi2, 4))
print("max diff:", round(np.max(np.abs(pi1 - pi2)), 6))   # 两种算法应一致

def hitting_time(start, target, maxsteps=20000, seed=0):
    rng = np.random.default_rng(seed)
    s, t = start, 0
    while s != target and t < maxsteps:
        s = rng.choice(len(P), p=P[s]); t += 1
    return t

ht = np.array([hitting_time(0, 4, seed=i) for i in range(5000)])
print("E[首达 0->4] ~", round(ht.mean(), 1))

幂迭代与线性求解得到的平稳分布应逐位一致（差值在 1e-6 量级），首达时为正值。两种求法对照是验证马尔可夫链实现正确性的强约束：幂迭代逼近极限，线性方程直接解特征方程，二者吻合即确认链可约且平稳分布存在。

七、鞅：公平博弈与可选停时

鞅要求 E[X_{t+1} | F_t] = X_t——在已知历史下，未来"无偏"。最经典的例子是布朗运动的 M_t = W_t² - t，由伊藤微分 d(W_t²) = 2W_t dW_t + dt 立得它是鞅。鞅的"无套利"直觉贯穿金融定价、Q-learning 的 TD 误差与收敛分析。

import numpy as np

def brownian_martingale(n=10000, dt=1e-3, seed=0):
    rng = np.random.default_rng(seed)
    dW = rng.normal(0, np.sqrt(dt), size=n)
    W = np.cumsum(dW)
    return W ** 2 - np.arange(1, n + 1) * dt     # M_t = W_t^2 - t

# 单条路径的 M_t 不是 0；鞅的"期望恒定"需在多条独立路径上取期望
rng = np.random.default_rng(7)
vals = np.array([brownian_martingale(seed=int(rng.integers(0, 10**9)))
                 for _ in range(2000)])
for i, label in enumerate(["t=1", "t=5", "t=10"]):
    print(f"E[M_t] @ {label} =", round(vals[:, i].mean(), 3))   # 各时刻都应≈0

M_t 在 t=1/5/10 处的无条件期望都应逼近 0——鞅的"期望恒定"性质。顺带一提：名为"鞅"的翻倍赌注策略之所以在公平赌局里期望为 0 却仍会破产，正是因为可选停时定理要求停时有界或满足可积条件，而无限翻倍破坏了这个前提——鞅的优雅与危险都在于此。

八、高斯过程：无限维高斯先验与非参数回归

高斯过程把"函数"本身当作随机变量：f ~ GP(m(x), k(x, x'))，任意有限点的联合分布都是高斯。给定带噪观测，后验仍是高斯，可直接给出预测均值与不确定性——这是贝叶斯深度学习、主动学习与不确定性量化的核心先验。

import numpy as np

def rbf(X, Z, ls=1.0, var=1.0):
    Xn, Zn = X / ls, Z / ls
    sq = (np.sum(Xn ** 2, 1)[:, None] + np.sum(Zn ** 2, 1)[None, :] - 2 * Xn @ Zn.T)
    return var * np.exp(-0.5 * np.clip(sq, 0, None))

def gp_posterior(X, y, Xs, noise=1e-3, ls=1.0, var=1.0):
    K = rbf(X, X, ls, var) + noise * np.eye(len(X))
    Ks = rbf(X, Xs, ls, var)
    Kss = rbf(Xs, Xs, ls, var)
    L = np.linalg.cholesky(K)
    alpha = np.linalg.solve(L.T, np.linalg.solve(L, y))
    mean = Ks.T @ alpha
    v = np.linalg.solve(L, Ks)
    cov = Kss - v.T @ v
    return mean, np.sqrt(np.diag(cov))

rng = np.random.default_rng(0)
X = np.linspace(0, 6, 8)[:, None]
y = np.sin(X[:, 0]) + rng.normal(0, 0.1, size=8)
Xs = np.linspace(0, 6, 200)[:, None]
mean, sd = gp_posterior(X, y, Xs)
mean_tr, _ = gp_posterior(X, y, X, noise=1e-3)
print("max|pred(train)-y| =", round(np.max(np.abs(mean_tr - y)), 3))  # 应很小（插值）
print("置信带宽度 ~", round(2 * sd.mean(), 3))                          # 训练区窄、外推区宽

后验均值在训练点上应几乎复现 y（误差在噪声量级），而置信带在观测密集处窄、外推处宽——不确定性随距离自然增长，这正是 GP 相较点估计模型最独特的价值：它给出可解释的置信，而非只会吐一个值的黑箱。

九、与 AI / 深度学习的映射表

十、工程边界与陷阱

• 离散化误差：伊藤 SDE 的欧拉-丸山法强收敛 O(√dt)、弱收敛 O(dt)；想要精确样本路径需小步长或高阶格式（Milstein）。
• 蒙特卡洛方差：均值估计误差 O(1/√M)，要降 10 倍方差需 100 倍样本；优先用对偶变量、控制变量、分层抽样。
• 平稳性假设：非平稳数据流（概念漂移）下马尔可夫/平稳假设失效，必须在线检测或更换窗口。
• 维数灾难：高斯过程随点数 O(n³)、高维 SDE 几乎不可解；高维场景转向稀疏 GP、变分推断或神经参数化。
• 自相关与时间尺度：强自相关会大幅削减有效样本量（ESS），MCMC 链必须诊断混合与自相关时间。
• 可复现性：随机性依赖种子，生产环境务必集中管理与记录 RNG 状态，避免"跑不出论文数字"。

结论

随机过程不是概率论的装饰，而是把不确定性变成可计算、可验证、可部署结构的工程语言。从随机游走的 CLT、泊松的计数律、布朗的二次变差，到伊藤 SDE、马尔可夫平稳分布、鞅的无偏性与高斯过程的非参数后验——每一条主线都有封闭的数学刻画和纯 NumPy 的落地实现。理解它们，就等于拿到了 SGD、MCMC、扩散模型、PageRank 与贝叶斯不确定性这些 AI 基石的"源代码"。下一个真正值得深挖的角度，是把这些过程放到高维非线性系统里（随机偏微分方程、神经 SDE、非平衡稳态），那才是连接统计物理与现代生成式 AI 的最后一公里。