第 24b 章 马尔可夫链蒙特卡洛

本章接第 24a 章，对应 Wasserman 原书第 24 章的后半部分（24.4–24.5 节）。基本蒙特卡洛要求能从目标分布直接抽样，重要性抽样要求找到与目标相似的提议分布，两者在高维后验面前都会失效。**马尔可夫链蒙特卡洛（MCMC）**换了一个思路：不直接从目标分布 \(f\) 抽样，而是构造一条以 \(f\) 为平稳分布的马尔可夫链，让它跑足够长，再把链上的点当作来自 \(f\) 的（相关的）样本。第 23 章的遍历定理保证时间平均收敛到 \(f\) 下的期望，细致平衡条件告诉我们怎样构造这样的链。MCMC 让贝叶斯方法在随机波动率模型、分层模型、贝叶斯 VAR 等量化模型中变得可行。

学习目标

1. 理解 MCMC 的基本思想：用遍历马尔可夫链的时间平均代替独立样本平均。
2. 会写出 Metropolis–Hastings 算法，并用细致平衡证明 \(f\) 是其平稳分布；理解为什么只需知道 \(f\) 的比例形式。
3. 理解提议尺度对混合的影响，会用接受率、轨迹图、自相关与有效样本量诊断链的质量。
4. 掌握随机游走 MH、独立 MH、Gibbs 抽样和 Metropolis-within-Gibbs 的适用场景。
5. 会推导正态分层模型的全条件分布并实现 Gibbs 抽样，理解分层模型的"收缩"效应及其在 alpha 估计中的应用。

读前导读

这一章在解决什么问题。 贝叶斯方法的结论是一个后验分布：参数的所有可能取值以及各自的可信度。要用它，几乎总是要算后验下的期望，比如"alpha 的后验均值""alpha 大于 0 的后验概率"。这些都是积分。参数只有一两个时可以数值积分；有几十个参数（40 位经理各一个 alpha，再加总体均值和离散度）时，积分算不动，也没法直接从后验抽样。第 24a 章的蒙特卡洛要求能直接抽样，重要性抽样要求找到一个和后验很像的分布，在高维下都会失效。

MCMC 的办法是：设计一个随机"漫步"规则，让一个点在参数空间里走，走的规则保证它在每个区域停留的时间比例恰好等于后验概率。走得足够久之后，把它走过的位置记下来，就相当于一堆来自后验的样本。你在 CFA 里用历史模拟法算 VaR，是拿一串历史收益当样本；MCMC 是自己造一串"会收敛到目标分布"的样本。区别在于这串样本前后相关，所以 2 万个样本不等于 2 万份信息，这就是本章反复强调"有效样本量"的原因。

本章的量化落点很直接：40 位基金经理的样本 alpha 噪声很大，分层模型把每个人的估计往"大家的平均"拉，拉多少由各自数据的可信度决定。这和 Black–Litterman 把观点向均衡收益收缩、信用评级里把小样本违约率向行业均值修正，是同一种思想。

需要先想起来的数学。

• 条件概率与贝叶斯公式（比例形式）：后验 \(\propto\) 似然 \(\times\) 先验，写成 \(f(\theta\mid x)\propto\mathcal L(\theta)f(\theta)\)。符号 \(\propto\) 读作"正比于"，意思是只差一个不依赖 \(\theta\) 的常数。例如 \(f(\theta)\propto e^{-\theta^2/2}\) 就是标准正态，常数 \(1/\sqrt{2\pi}\) 被省掉了。本章的核心技巧就是"常数不知道也没关系"。参见 第 00 册第 07 章 概率中的分析工具。
• 期望就是积分：\(\mathbb E_f h(X)=\int h(x)f(x)dx\)，下标 \(f\) 表示"在 \(X\) 服从 \(f\) 时求期望"。取 \(h(x)=I(|x|<1)\)（示性函数，条件成立取 1 否则取 0），期望就变成概率 \(\mathbb P(|X|<1)\)。所以"估计概率"和"估计期望"是同一件事。参见 第 00 册第 03 章 积分。
• 马尔可夫链与平稳分布（第 23 章）：马尔可夫链是"下一步只看当前位置"的随机过程。平稳分布是一种状态分布，链按规则走一步之后分布不变。类比信用评级迁移矩阵：如果某个评级分布乘以迁移矩阵后还是自己，它就是平稳分布。遍历定理说：链走得够久，时间平均等于平稳分布下的期望。
• 正态密度的"配方"：正态密度取对数后是 \(x\) 的二次函数 \(-\frac{(x-m)^2}{2v}\)。反过来，若某个密度取对数后是 \(-\frac{a}{2}x^2+bx+\text{常数}\)，它就是均值 \(b/a\)、方差 \(1/a\) 的正态。例：\(\log f=-x^2+4x\)，则 \(a=2,b=4\)，是 \(N(2,0.5)\)。这是推导全条件分布的主要工具。
• 变量替换与 Jacobian（一维）：若 \(W=\log X\)，则 \(f_W(w)=f_X(e^w)\cdot\big|\frac{dx}{dw}\big|=f_X(e^w)e^w\)。多出的因子 \(e^w\) 是导数，用来修正"刻度被拉伸"。参见 第 00 册第 05 章 多元微积分与优化 中雅可比一节，以及 第 02 章 导数与泰勒展开。

怎么读这一章。 24b.1.1–24b.1.3 是核心必读：算法本身和它为什么对。24b.1.4–24b.1.5 讲实际用起来要注意什么，务必读，尤其是有效样本量。24b.2.1–24b.2.3 介绍几种变体，先抓住"随机游走 MH 调步长""Gibbs 一次只动一个坐标"两个要点即可。24b.2.4 的分层模型推导较密，第一次可以只看结论（精度加权平均与收缩），读完 24b.3 的基金经理例子后再回头细看推导。24b.2.5 第一次可略读。24b.3 的解读第 2、3 条是本章对量化工作最有价值的部分。

24b.1 Metropolis–Hastings 算法

24b.1.1 思想

我们要计算 \(I=\int h(x)f(x)dx\)，但无法从 \(f\) 直接抽样。若能构造一条马尔可夫链 \(X_1,X_2,\dots\)，其平稳分布恰为 \(f\)，则在一定条件下（第 23 章定理 23.25）

白话解释：\(\xrightarrow{P}\) 是"依概率收敛"：\(N\) 越大，左边的平均值偏离 \(I\) 超过任意给定小量的概率越接近 0。这和大数定律是同一种结论，只是大数定律要求样本独立，这里样本来自一条链、前后相关，所以需要第 23 章的遍历定理来替代。"一定条件"主要指链不会被困在某个局部（不可约）、不会周期性地来回跳（非周期）等，MH 算法在常见设置下都满足。

24b.1.2 算法

选一个提议分布（proposal distribution） \(q(y\mid x)\)：任何易于抽样的条件密度。

Metropolis–Hastings 算法。 任选初值 \(X_0\)。已得到 \(X_0,\dots,X_i\) 后：

1. 生成候选 \(Y\sim q(y\mid X_i)\)；
2. 计算接受概率 \(r=r(X_i,Y)\)，
   \[r(x,y)=\min\Big\{\frac{f(y)}{f(x)}\frac{q(x\mid y)}{q(y\mid x)},\ 1\Big\};\]
3. 以概率 \(r\) 令 \(X_{i+1}=Y\)（接受），以概率 \(1-r\) 令 \(X_{i+1}=X_i\)（拒绝，原地不动）。

实现第 3 步：抽 \(U\sim\text{Uniform}(0,1)\)，若 \(U<r\) 则接受（注 24.8）。常用的提议是 \(q(y\mid x)=N(x,b^2)\)，它关于 \(x,y\) 对称，此时 \(r=\min\{f(Y)/f(X_i),1\}\)（注 24.9）：往概率更高的地方走总是接受，往概率更低的地方走以概率 \(f(y)/f(x)\) 接受。

关键性质：\(f\) 只出现在比值 \(f(y)/f(x)\) 中，所以只需知道 \(f\) 到一个比例常数。 对贝叶斯后验 \(f(\theta\mid x)\propto\mathcal L(\theta)f(\theta)\)，那个最难算的归一化常数 \(c=\int\mathcal L f\) 根本不需要。实现时为避免数值下溢，应在对数尺度上比较：\(\log U<\log f(Y)-\log f(X_i)\)。

白话解释：为什么归一化常数会消掉？若 \(f(\theta)=\mathcal L(\theta)f_0(\theta)/c\)（\(f_0\) 为先验），那么 \(\frac{f(y)}{f(x)}=\frac{\mathcal L(y)f_0(y)/c}{\mathcal L(x)f_0(x)/c}=\frac{\mathcal L(y)f_0(y)}{\mathcal L(x)f_0(x)}\)，\(c\) 上下约掉。贝叶斯计算里最难的就是 \(c=\int\mathcal Lf_0\,d\theta\) 这个高维积分，MH 完全绕开了它。这和比较两只股票的相对估值类似：算市盈率之比时，不需要知道整个市场的总市值。

小数值例：目标为标准正态，只用 \(\log f(x)=-x^2/2\)。当前 \(x=0.5\)，候选 \(y=1.5\)，则 \(\log f(y)-\log f(x)=-1.125+0.125=-1\)，接受概率 \(r=e^{-1}\approx0.37\)；若候选 \(y=0.2\)，差值为正，直接接受。对数尺度的理由是：几十个观测的似然连乘后可能是 \(10^{-300}\) 这种量级，直接相除会下溢成 \(0/0\)，取对数后只是两个普通数相减。

24b.1.3 为什么有效：细致平衡

用 \(p(x,y)\) 表示从 \(x\) 转移到 \(y\) 的转移密度。\(f\) 是平稳分布意味着 \(f(x)=\int f(y)p(y,x)dy\)。与第 23 章离散情形相同，细致平衡

证明 MH 满足细致平衡。 取 \(x\ne y\)，不妨设 \(f(x)q(y\mid x)>f(y)q(x\mid y)\)（相等的情形概率为零）。则

推导拆解：分三步看。

第一步，细致平衡为什么推出平稳。平稳的定义是"走一步后分布不变"：到达 \(x\) 的总概率 \(\int f(y)p(y,x)dy\) 等于 \(f(x)\)。细致平衡把被积函数 \(f(y)p(y,x)\) 换成 \(f(x)p(x,y)\)；\(f(x)\) 与积分变量 \(y\) 无关，提到积分外；剩下 \(\int p(x,y)dy\) 是"从 \(x\) 出发走到某处"的总概率，等于 1。细致平衡是比平稳更强的条件：它要求每一对 \(x,y\) 之间的往返流量都相等，而平稳只要求每个点的总流入等于总流出。

第二步，"不妨设"。\(f(x)q(y\mid x)\) 与 \(f(y)q(x\mid y)\) 必有一个大，另一种情况把 \(x,y\) 对调即可，证明完全对称，所以只证一种。由假设，MH 比值 \(\frac{f(y)q(x\mid y)}{f(x)q(y\mid x)}<1\)，取 min 后就是它本身；而反方向的比值是它的倒数，大于 1，取 min 后为 1。

第三步，转移密度 \(p(x,y)\) 为什么是 \(q(y\mid x)r(x,y)\)。从 \(x\) 移到另一个点 \(y\neq x\)，必须"先提议到 \(y\)"（密度 \(q(y\mid x)\)）且"被接受"（概率 \(r\)），两件事相乘。拒绝时停在原地，这部分概率只落在 \(y=x\) 上，所以证明里只考虑 \(x\ne y\) 即可，\(x=y\) 时细致平衡两边相同，自动成立。

金融直觉：细致平衡像两个账户之间的资金往来轧差为零：\(x\) 账户流向 \(y\) 的金额 \(f(x)p(x,y)\) 等于 \(y\) 流回 \(x\) 的金额。每一对账户都轧平，所有账户的余额（分布）自然不变。

接受概率中的 \(\frac{q(x\mid y)}{q(y\mid x)}\) 正是为了抵消提议分布本身的不对称：若提议更容易从 \(x\) 走到 \(y\) 而不容易走回来，就要相应降低接受概率，使两个方向的概率流相等。

24b.1.4 例 24.10：提议尺度与混合

目标为 Cauchy 密度 \(f(x)=\frac1\pi\frac1{1+x^2}\)，提议 \(N(x,b^2)\)，\(r(x,y)=\min\{\frac{1+x^2}{1+y^2},1\}\)。原书画了三条长 1000 的链：

• \(b=0.1\)：步子太小，几乎每步都被接受，但链移动缓慢，像在原地徘徊，长时间探索不到尾部，直方图与真密度差距大；
• \(b=10\)：步子太大，提议经常落到远处的尾部，\(r\) 很小，频繁被拒绝，链长时间"卡"在同一个值；
• \(b=1\)：居中，链更快地覆盖目标分布。

若链的样本很快"看起来像"目标分布，称链混合良好（mixing well）。构造混合良好的链有一定艺术性。

24b.1.5 实践诊断（本教材补充）

原书点到为止，实际使用 MCMC 时以下几项是必做的：

• 预烧期（burn-in）：初值可能远离分布的主体，丢弃前若干次迭代。
• 轨迹图（trace plot）：把 \(X_i\) 对 \(i\) 作图。混合良好的链看起来像"毛毛虫"，在稳定范围内快速上下波动；混合差的链有长时间的平台（频繁拒绝）或缓慢漂移（步子太小）。
• 自相关与有效样本量（ESS）：若链的自相关为 \(\rho_k\)，均值估计的方差约为 \(\frac{\sigma^2}N(1+2\sum_k\rho_k)\)，相当于 \(\text{ESS}=N/(1+2\sum_k\rho_k)\) 个独立样本。报告 MCMC 结果时应报告 ESS 而非迭代次数。

推导拆解：ESS 公式从哪来。设链已平稳，\(\mathbb V(X_i)=\sigma^2\)，\(\text{Cov}(X_i,X_{i+k})=\rho_k\sigma^2\)（\(\rho_k\) 为相隔 \(k\) 步的自相关）。样本均值的方差是所有两两协方差之和除以 \(N^2\)：\(\mathbb V(\bar X)=\frac1{N^2}\sum_i\sum_j\text{Cov}(X_i,X_j)\)。对角线 \(N\) 项各贡献 \(\sigma^2\)；相隔 \(k\) 的项有 \(N-k\) 对、上下两侧各一次。当 \(N\) 远大于自相关衰减所需的步数时，\(N-k\approx N\)，得到 \(\mathbb V(\bar X)\approx\frac{\sigma^2}N(1+2\sum_{k\ge1}\rho_k)\)。把它写成 \(\sigma^2/\text{ESS}\) 就得到 ESS 公式。这就是你熟悉的"组合方差 = 各项方差 + 2 倍协方差之和"在时间序列上的版本。

例：\(\rho_k=0.5^k\)，则 \(\sum_{k\ge1}0.5^k=1\)（几何级数），\(1+2\times1=3\)，1 万次迭代只值约 3,333 个独立样本。

24b.2 MCMC 的几种变体

24b.2.1 随机游走 MH

提议 \(Y=X_i+\epsilon_i\)，\(\epsilon_i\sim g\)，\(g\) 对称（常取 \(N(0,b^2)\)），于是 \(q(y\mid x)=g(y-x)\)，\(r=\min\{1,f(y)/f(x)\}\)。如果没有接受–拒绝这一步，链就是普通的随机游走，故名。难点是选 \(b\)。

经验法则：选 \(b\) 使提议约有 50% 被接受。（补充：Roberts、Gelman 等人的理论结果表明，在高维近似正态的目标下最优接受率约为 23%，一维时约为 44%。）

警告：随机游走 MH 只在 \(X\) 取值于整个实轴时自然。若 \(X\) 受限（如方差 \(\sigma^2>0\)），应先变换（如 \(W=\log X\)），对 \(W\) 的分布做 MH，再变换回来。注意 \(W\) 的密度要乘上 Jacobian：若 \(X\) 的密度为 \(f_X\)，则 \(f_W(w)=f_X(e^w)e^w\)。

推导拆解：Jacobian 因子的来源。\(W\le w\) 等价于 \(X\le e^w\)，所以 \(F_W(w)=F_X(e^w)\)。两边对 \(w\) 求导，右边用链式法则（外层 \(F_X\) 的导数是 \(f_X\)，内层 \(e^w\) 的导数是 \(e^w\)），得 \(f_W(w)=f_X(e^w)\cdot e^w\)。在对数尺度上实现时就是 \(\log f_W(w)=\log f_X(e^w)+w\)，多加一个 \(w\)。

漏掉这一项会怎样：链会对 \(W\) 抽样，但目标错了，变换回 \(X\) 后得到的分布不是 \(f_X\)。如果直接对 \(\sigma\) 加扰动，提议会落到 \(\sigma\le0\)；此时 \(f=0\)、必被拒绝，算法仍然正确，但在 0 附近浪费大量提议，而且代码里 \(\log\sigma\) 会报错，这是实践中更常见的问题。

24b.2.2 独立 MH

提议来自一个与当前状态无关的固定分布 \(g\)（通常是 \(f\) 的近似，比如以后验众数为中心、以 Fisher 信息的逆为协方差的 \(t\) 分布）：

24b.2.3 Gibbs 抽样

前两种方法原则上适用于任意维数，但高维时很难调到混合良好。Gibbs 抽样把高维问题拆成一系列一维（或低维）问题。设 \((X,Y)\sim f_{X,Y}\)，假设能从两个全条件分布 \(f_{X\mid Y}\)、\(f_{Y\mid X}\) 直接抽样。从 \((X_0,Y_0)\) 出发，

（补充：Gibbs 的每一步都可以看作提议分布为全条件分布的 MH 步，代入 MH 接受概率公式可以验证 \(r\equiv1\)——提议总是被接受。这就是 Gibbs 正确性的来源。）

推导拆解：以二维为例，当前点 \((x,y)\)，Gibbs 提议只改 \(x\)：候选 \((x',y)\)，提议密度 \(q(x'\mid x)=f(x'\mid y)\)，反方向 \(q(x\mid x')=f(x\mid y)\)。代入 MH 比值：

\[\frac{f(x',y)}{f(x,y)}\cdot\frac{f(x\mid y)}{f(x'\mid y)}=\frac{f(x'\mid y)f(y)}{f(x\mid y)f(y)}\cdot\frac{f(x\mid y)}{f(x'\mid y)}=1.\]

第一个等号用了"联合 = 条件 × 边际"：\(f(x,y)=f(x\mid y)f(y)\)。所有因子两两约掉，接受概率为 \(\min\{1,1\}=1\)。直观上，全条件分布本身就是"在 \(y\) 固定时 \(x\) 的正确分布"，从它抽出来的点不需要再纠偏。

白话解释："全条件分布"\(f(x_j\mid x_{-j})\) 中的 \(x_{-j}\) 是记号习惯，表示"除第 \(j\) 个以外的所有坐标"，下文的 \(\text{rest}\) 是同一个意思。"共轭"指先验和似然属于配套的分布族，乘起来后验仍在同一族（如正态先验配正态似然得正态后验），于是全条件分布有现成的抽样函数。

Gibbs 的前提是全条件分布容易抽样，这在共轭结构下往往成立。它的弱点是：若参数之间高度相关，逐坐标移动会很慢（在狭长的山脊上"之"字形前进）。

24b.2.4 例 24.11：正态分层模型

问题。 抽 \(k\) 个城市，第 \(i\) 个城市调查 \(n_i\) 人，\(Y_i\) 人患病，\(Y_i\sim\text{Binomial}(n_i,p_i)\)。各城市患病率不同，把 \(p_1,\dots,p_k\) 看作来自某个总体分布 \(F\) 的随机抽样。目标是估计每个 \(p_i\) 和总体平均患病率。

正态近似。 对 logit 做 Delta 方法：令 \(\psi_i=\log\frac{p_i}{1-p_i}\)，\(Z_i=\hat\psi_i=\log\frac{\hat p_i}{1-\hat p_i}\)，则 \(Z_i\approx N(\psi_i,\sigma_i^2)\)，\(\sigma_i^2=\frac1{n_i\hat p_i(1-\hat p_i)}\)（logit 尺度上的正态近似比原尺度更准），并把 \(\sigma_i\) 当作已知。

白话解释：这一步只是在准备数据。比例 \(\hat p_i\) 限制在 0 到 1 之间，且方差随 \(p\) 变化，不便套正态模型；logit 变换 \(\log\frac p{1-p}\)（即"对数赔率"）把它拉到整个实轴上。Delta 方法是一阶泰勒近似：\(g(\hat p)\approx g(p)+g'(p)(\hat p-p)\)，所以 \(\mathbb V(g(\hat p))\approx g'(p)^2\mathbb V(\hat p)\)。这里 \(g'(p)=\frac1{p(1-p)}\)，\(\mathbb V(\hat p)=\frac{p(1-p)}n\)，相乘得 \(\frac1{np(1-p)}\)。这和用久期近似债券价格变动是同一个思路：价格变化 ≈ 导数 × 收益率变化。

模型为

全条件分布。 只保留含目标参数的因子，配方：

• \(f(\mu\mid\text{rest})\propto\prod_i\exp\{-\frac12(\psi_i-\mu)^2\}\propto\exp\{-\frac k2(\mu-b)^2\}\)，\(b=\bar\psi\)，即 \(\mu\mid\text{rest}\sim N(\bar\psi,1/k)\)；
• \(f(\psi_i\mid\text{rest})\propto\exp\{-\frac12(\psi_i-\mu)^2\}\exp\{-\frac1{2\sigma_i^2}(Z_i-\psi_i)^2\}\)，两个正态核相乘仍为正态：
  \[\psi_i\mid\text{rest}\sim N(e_i,d_i^2),\qquad e_i=\frac{Z_i/\sigma_i^2+\mu}{1/\sigma_i^2+1},\qquad d_i^2=\frac1{1/\sigma_i^2+1}.\]

\(e_i\) 是数据 \(Z_i\) 与总体均值 \(\mu\) 的精度加权平均：样本小（\(\sigma_i^2\) 大）的城市被更多地拉向总体均值，样本大的城市基本保持自己的数据。

推导拆解：两个全条件分布的配方过程。

\(\mu\) 的全条件：后验里只有第一组因子含 \(\mu\)，第二组 \(\exp\{-\frac1{2\sigma_i^2}(Z_i-\psi_i)^2\}\) 与 \(\mu\) 无关，当作常数丢掉。指数部分 \(-\frac12\sum_i(\psi_i-\mu)^2\) 展开为 \(-\frac k2\mu^2+\mu\sum_i\psi_i+\text{常数}\)。按导读里的规则，\(a=k\)、\(b=\sum_i\psi_i\)，均值 \(b/a=\bar\psi\)，方差 \(1/a=1/k\)。直观上就是"把 \(k\) 个 \(\psi_i\) 当作来自 \(N(\mu,1)\) 的样本估计 \(\mu\)"，平坦先验下后验以样本均值为中心、方差为 \(1/k\)。

\(\psi_i\) 的全条件：只保留含 \(\psi_i\) 的两个因子，指数部分为 \(-\frac12(\psi_i-\mu)^2-\frac1{2\sigma_i^2}(\psi_i-Z_i)^2\)。\(\psi_i^2\) 的系数为 \(-\frac12(1+\frac1{\sigma_i^2})\)，故 \(a=1+1/\sigma_i^2\)；\(\psi_i\) 的一次项系数为 \(\mu+Z_i/\sigma_i^2=b\)。均值 \(b/a\) 即 \(e_i\)，方差 \(1/a\) 即 \(d_i^2\)。一般 \(\tau\) 时把"1"换成 \(1/\tau^2\)，就是本章小结表中的公式。

数值例：\(Z_i=2\)，\(\sigma_i^2=1\)，\(\mu=0\)，\(\tau=1\)，则 \(e_i=(2+0)/(1+1)=1\)，向均值拉回一半；若 \(\sigma_i^2=0.25\)（样本大 4 倍），\(e_i=(8+0)/(4+1)=1.6\)，只拉回 20%。

金融直觉：这就是"精度加权"，精度即方差的倒数。它和最小方差组合给每项资产按 \(1/\sigma^2\) 分配权重是同一个原理：两个对同一量的独立估计，按各自的精度加权合并，合并后方差最小，且合并后的精度等于两个精度之和（\(1/d_i^2=1/\sigma_i^2+1/\tau^2\)）。分析师给小盘股定盈利预测时，公司自身历史数据短就更多参考行业均值，也是同一做法。

Gibbs 迭代：抽 \(\mu\sim N(\bar\psi,1/k)\)；依次抽 \(\psi_i\sim N(e_i,d_i^2)\)，\(i=1,\dots,k\)；重复 \(N\) 次。最后用 \(p_i=e^{\psi_i}/(1+e^{\psi_i})\) 变换回概率尺度。

原书用 \(k=20\) 个城市、每城 20 人的数据演示：\(p_1\) 与 \(\mu\) 的轨迹图显示链混合良好；Bayes 估计相对于原始比例被**收缩（shrunk）**到一起。\(\tau\) 控制收缩程度——\(\tau\) 越小，各城市越被拉向共同均值。原书固定 \(\tau=1\)，并指出实践中应把 \(\tau\) 当作未知参数由数据决定（本章量化实战会这样做）。

24b.2.5 Metropolis-within-Gibbs

若某个全条件分布无法直接抽样，就对该坐标做一步 MH。二维情形，\(q\)、\(\tilde q\) 分别为 \(x\)、\(y\) 的提议分布：

1. 抽 \(Z\sim q(z\mid X_n)\)，\(r=\min\Big\{\frac{f(Z,Y_n)}{f(X_n,Y_n)}\frac{q(X_n\mid Z)}{q(Z\mid X_n)},1\Big\}\)，以概率 \(r\) 令 \(X_{n+1}=Z\)，否则 \(X_{n+1}=X_n\)；
2. 抽 \(Z\sim\tilde q(z\mid Y_n)\)，\(r=\min\Big\{\frac{f(X_{n+1},Z)}{f(X_{n+1},Y_n)}\frac{\tilde q(Y_n\mid Z)}{\tilde q(Z\mid Y_n)},1\Big\}\)，以概率 \(r\) 令 \(Y_{n+1}=Z\)，否则 \(Y_{n+1}=Y_n\)。

对 \(X\) 做 MH 步时 \(Y\) 固定，反之亦然；可以推广到任意维，对能直接抽样的坐标用 Gibbs、不能的用 MH。贝叶斯 logistic 回归（原书习题 5）、随机波动率模型中的波动率路径，都常用这种混合方案。

24b.3 量化实战

场景。 (1) 重做例 24.10：随机游走 MH 抽 Cauchy 分布，比较不同提议尺度 \(b\) 的接受率、估计精度和有效样本量。(2) 例 24.11 的量化版本：40 位基金经理各有 60 个月的业绩，样本 alpha（月度，%）的抽样标准误由各自的残差波动决定。真实 alpha 来自 \(N(0.20,0.15^2)\)。用分层模型 \(\alpha_i\sim N(\mu,\tau^2)\)、\(Z_i\mid\alpha_i\sim N(\alpha_i,s_i^2)\) 做 Gibbs 抽样，并且把 \(\tau\) 也当作未知参数（取 \(\tau\) 的均匀先验，\(\tau^2\) 的全条件分布为逆 Gamma）。比较原始样本 alpha 与收缩估计的误差，并考察样本 alpha 最高的 5 位经理。

推导拆解：\(\tau^2\) 的全条件为什么是逆 Gamma（IG，定义见附录 A1.7）。含 \(\tau\) 的因子只有 \(\prod_i\frac1\tau\exp\{-\frac{(\psi_i-\mu)^2}{2\tau^2}\}=(\tau^2)^{-k/2}\exp\{-\frac S{2\tau^2}\}\)，其中 \(S=\sum_i(\psi_i-\mu)^2\)。注意正态密度前面的 \(1/\tau\) 这次不能丢，因为它含 \(\tau\)。"\(\tau\) 均匀"换成 \(\tau^2\) 的先验时要乘 Jacobian：\(\tau=\sqrt{\tau^2}\)，\(\frac{d\tau}{d\tau^2}=\frac1{2\tau}\)，所以 \(\tau^2\) 的先验 \(\propto(\tau^2)^{-1/2}\)。两者相乘得 \((\tau^2)^{-(k+1)/2}\exp\{-\frac{S/2}{\tau^2}\}\)，对照 IG\((a,b)\) 的核 \(x^{-a-1}e^{-b/x}\)，得 \(a=\frac{k-1}2\)、\(b=\frac S2\)。代码中 (S/2)/rng.gamma((k-1)/2) 正是利用"IG 变量 = \(b\) 除以 Gamma\((a,1)\) 变量"来抽样。

import numpy as np
rng = np.random.default_rng(2424)

def ess(x):                                          # 有效样本量：N / (1 + 2Σρ_k)，累加到自相关首次为负
    x = x - x.mean(); n = len(x)
    f = np.fft.rfft(x, 2 * n); ac = np.fft.irfft(f * np.conj(f))[:n]; ac /= ac[0]
    s = 0.0
    for k in range(1, n):
        if ac[k] < 0: break
        s += ac[k]
    return n / (1 + 2 * s)

# ---------- (1) 例 24.10：随机游走 MH 抽 Cauchy，提议尺度 b 的影响 ----------
log_f = lambda x: -np.log1p(x ** 2)                  # Cauchy 密度，只需到比例常数
N = 20000
print("   b    接受率   P(|X|<1) 估计(真值0.5)   有效样本量ESS")
for b in [0.1, 1.0, 2.5, 10.0]:
    x = np.empty(N); x[0] = 0.0; acc = 0
    for i in range(N - 1):
        y = x[i] + b * rng.standard_normal()
        if np.log(rng.random()) < log_f(y) - log_f(x[i]):
            x[i + 1] = y; acc += 1
        else:
            x[i + 1] = x[i]
    ind = (np.abs(x) < 1).astype(float)
    print("%5.1f  %6.2f   %8.3f            %7.0f" % (b, acc / (N - 1), ind.mean(), ess(ind)))

# ---------- (2) 例 24.11 的量化版本：基金经理 alpha 的分层收缩（Gibbs，τ 也作为未知参数） ----------
k, T = 40, 60                                        # 40 位经理，各 60 个月
mu_true, tau_true = 0.20, 0.15                       # 月度 alpha（%）的总体均值与离散度
alpha = rng.normal(mu_true, tau_true, k)
sig = rng.uniform(1.0, 4.0, k)                       # 各经理的月度残差波动（%）
Z = alpha + sig / np.sqrt(T) * rng.standard_normal(k) # 观测到的样本 alpha
s2 = (sig / np.sqrt(T)) ** 2                          # 视为已知的抽样方差 σ_i²
t_stat = Z / np.sqrt(s2)

iters, burn = 6000, 1000
mu, tau2, psi = Z.mean(), Z.var(), Z.copy()
keep_mu, keep_tau, keep_psi = [], [], []
for it in range(iters):
    prec = 1 / tau2 + 1 / s2                          # ψ_i | rest ~ N(e_i, d_i²)，精度加权
    psi = rng.normal((Z / s2 + mu / tau2) / prec, np.sqrt(1 / prec))
    mu = rng.normal(psi.mean(), np.sqrt(tau2 / k))    # μ | rest ~ N(ψ̄, τ²/k)（平坦先验）
    S = np.sum((psi - mu) ** 2)                       # τ 取均匀先验 ⇒ τ² | rest ~ IG((k-1)/2, S/2)
    tau2 = (S / 2) / rng.gamma((k - 1) / 2)
    if it >= burn:
        keep_mu.append(mu); keep_tau.append(np.sqrt(tau2)); keep_psi.append(psi)
keep_psi = np.array(keep_psi); post = keep_psi.mean(0)
print("μ 后验均值=%.3f (真值 %.2f)  95%%区间=(%.3f, %.3f)" % (np.mean(keep_mu), mu_true, *np.quantile(keep_mu, [0.025, 0.975])))
print("τ 后验均值=%.3f (真值 %.2f)  95%%区间=(%.3f, %.3f)" % (np.mean(keep_tau), tau_true, *np.quantile(keep_tau, [0.025, 0.975])))
print("均方误差 (×1e3)：原始样本alpha=%.2f   分层收缩后=%.2f" % (1e3 * np.mean((Z - alpha) ** 2), 1e3 * np.mean((post - alpha) ** 2)))
top = np.argsort(Z)[-5:][::-1]
print("样本 alpha 最高的 5 位：")
print("  样本alpha  t值   收缩后  P(α>0)  真实alpha")
for i in top:
    print("  %7.2f  %5.1f  %6.2f  %6.2f  %7.2f" % (Z[i], t_stat[i], post[i], (keep_psi[:, i] > 0).mean(), alpha[i]))

关键输出：

   b    接受率   P(|X|<1) 估计(真值0.5)   有效样本量ESS
  0.1    0.97      0.661                 41
  1.0    0.77      0.525                578
  2.5    0.56      0.509               2061
 10.0    0.27      0.512               1497
μ 后验均值=0.187 (真值 0.20)  95%区间=(0.081, 0.292)
τ 后验均值=0.160 (真值 0.15)  95%区间=(0.023, 0.315)
均方误差 (×1e3)：原始样本alpha=99.50   分层收缩后=14.94
样本 alpha 最高的 5 位：
  样本alpha  t值   收缩后  P(α>0)  真实alpha
     0.91    1.9    0.27    0.95     0.24
     0.89    4.2    0.43    1.00     0.27
     0.79    1.6    0.25    0.95     0.18
     0.76    1.5    0.24    0.94     0.26
     0.74    2.0    0.28    0.97     0.49

解读。

1. 接受率不是越高越好。 \(b=0.1\) 时接受率 97%，但 2 万次迭代只相当于 41 个独立样本，\(\mathbb P(|X|<1)\) 的估计 0.661 偏得离谱——链根本没走到 Cauchy 的厚尾。接受率 56%（\(b=2.5\)）时 ESS 最大，约为迭代次数的 10%，与"约 50% 接受率"的经验法则一致。\(b=10\) 接受率只有 27%，ESS 也下降。报告 MCMC 结果时，迭代次数不说明任何问题，ESS 才说明问题。
2. 收缩大幅降低估计误差。 分层模型把 40 个样本 alpha 的均方误差从 0.0995 降到 0.0149，约为原来的 15%。原因是样本 alpha 的抽样噪声（标准误 0.13%–0.52%）远大于经理之间真实 alpha 的离散度（0.15%），原始估计主要反映的是噪声。
3. "冠军"往往是运气。 样本 alpha 最高的经理（月度 0.91%，年化约 11%）真实 alpha 只有 0.24%；他的 t 值只有 1.9，模型把他收缩到 0.27%。t 值最高（4.2）的第二名波动小、数据可信，收缩得较少（0.43%），但仍被大幅下调。真实 alpha 最高（0.49%）的经理反而排在第五。若按样本 alpha 排序配置资金，选中的多半是高波动、靠运气领先的经理。这与第 12 章 12.8 节的"赢家诅咒"、第 10b 章的数据挖掘偏差是同一现象，分层贝叶斯给出了系统的修正方法；Black–Litterman 模型中把观点向均衡收益收缩，也是同样的思想（第 11 册）。
4. \(\tau\) 的不确定性很大。 \(\tau\) 的 95% 后验区间为 \((0.023,0.315)\)：40 个经理不足以精确估计"真实 alpha 有多分散"。固定 \(\tau\)（像原书那样）会低估这部分不确定性；把 \(\tau\) 当作参数让后验自动反映它，这正是全贝叶斯处理的优点。

金融直觉：解读第 3 条中的 \(P(\alpha>0)\) 一列可以和频率派的 t 检验对照着看。t 值 1.9 的经理在 5% 双侧检验下不显著，但后验 \(P(\alpha>0)=0.95\) 看起来很"确定"。两者回答的问题不同：p 值说的是"若真实 alpha 为 0，看到这么高样本 alpha 的概率"；后验概率说的是"在整个经理群体的先验下，这位经理 alpha 为正的可信度"。这里群体均值本身就是正的（约 0.19%），所以几乎所有经理的 \(P(\alpha>0)\) 都会偏高，真正区分经理的是收缩后的点估计，而不是这一列。

MCMC 在量化中的其他应用。 随机波动率模型（潜在的对数波动率路径用 Gibbs 或粒子方法逐段更新）；贝叶斯 VAR 与因子模型；带参数不确定性的组合优化（对后验样本逐个求最优权重，再平均或取稳健解）；以及第 24a 章所说的，任何参数复杂函数（如 Sharpe 比、最大回撤的期望）的后验分布。

本章小结

MCMC 构造一条以目标分布 \(f\) 为平稳分布的遍历马尔可夫链，用时间平均估计 \(f\) 下的期望。Metropolis–Hastings 以概率 \(\min\{\frac{f(y)q(x\mid y)}{f(x)q(y\mid x)},1\}\) 接受提议，由细致平衡保证 \(f\) 平稳，且只需知道 \(f\) 的比例形式。随机游走 MH 的效率取决于提议尺度：步子太小链移动慢，步子太大频繁拒绝，接受率约 50%（高维约 23%）时较好；受限参数要先变换并带上 Jacobian。独立 MH 相当于 MCMC 版的重要性抽样。Gibbs 抽样逐坐标从全条件分布抽样，接受率恒为 1，适合共轭结构；不能直接抽样的坐标用 Metropolis-within-Gibbs。分层模型的 Gibbs 抽样产生精度加权的收缩估计，在噪声大、个体多的问题（基金 alpha、因子溢价、各资产预期收益）中能大幅降低估计误差。使用 MCMC 必须做诊断：预烧期、轨迹图、自相关、有效样本量、多链比较。

练习

基础

1. 证明 MH 算法满足细致平衡（不看书重写 24b.1.3 节的证明），并说明为什么不需要 \(f\) 的归一化常数。
2. 证明 Gibbs 抽样的每一步是接受率恒为 1 的 MH 步。（提示：提议 \(q(x_j'\mid x)=f(x_j'\mid x_{-j})\)，代入 MH 比值，利用 \(f(x)=f(x_j\mid x_{-j})f(x_{-j})\)。）
3. 对正态分层模型（一般 \(\tau\)），推导 \(\psi_i\) 和 \(\mu\) 的全条件分布；再在 \(\tau\) 的均匀先验下推导 \(\tau^2\) 的全条件分布。
4. 若对 \(\sigma>0\) 做随机游走 MH，直接对 \(\sigma\) 加正态扰动会出现什么问题？改为对 \(W=\log\sigma\) 做 MH 时，目标密度应如何修改？
5. 一条 MCMC 链的一阶自相关为 0.9，且自相关近似按 \(\rho_k=0.9^k\) 衰减。10 万次迭代的有效样本量约为多少？（答案：\(1+2\sum_k0.9^k=1+2\times9=19\)，ESS 约 5263。）

进阶

6. 原书习题 4：逆高斯分布 \(f(z)\propto z^{-3/2}\exp\{-\theta_1z-\theta_2/z+2\sqrt{\theta_1\theta_2}+\log\sqrt{2\theta_2}\}\)，\(\theta_1=1.5\)，\(\theta_2=2\)。(a) 用 Gamma 提议的独立 MH 抽 1000 个样本，与理论均值 \(\mathbb EZ=\sqrt{\theta_2/\theta_1}\)、\(\mathbb E(1/Z)=\sqrt{\theta_1/\theta_2}+\frac1{2\theta_2}\) 比较；(b) 令 \(W=\log Z\)，写出 \(W\) 的密度（含 Jacobian），用随机游走 MH 抽样后取指数。
7. 原书习题 5：对一个二分类数据集（可以用第 22a 章 22a.10 节的模拟涨跌数据：先运行那段代码，用其中的 make(n) 生成特征 Xtr 与标签 ytr）做贝叶斯 logistic 回归，平坦先验，用 Metropolis-within-Gibbs 抽 10,000 个后验样本，画各 \(\beta_j\) 的后验直方图，给出后验均值与 95% 后验区间，并与 MLE 和 Wald 区间比较。
8. 在 24b.3 节实验 (2) 中把经理人数从 40 改为 10 和 200，观察 \(\tau\) 的后验区间和收缩后均方误差如何变化。经理人数少时，固定 \(\tau\) 与估计 \(\tau\) 哪个更稳妥？
9. 用 Gibbs 抽样实现一个两状态马尔可夫切换的收益模型：给定状态路径，各状态的均值、方差有共轭后验；给定参数，用前向滤波–后向抽样（FFBS）抽取状态路径。这是第 23 章 regime 模型的贝叶斯估计方法。（综合题，可参考第 06 册。）

原书推荐习题：第 24 章习题 4、5。

原书对照