很多人第一次接触 MCMC，多半是被同一个问题逼到墙角的：手里有个分布，密度函数能算——甚至只有没归一化的密度核——但你就是没办法直接从里面“抽”出样本来。我当年做贝叶斯模型评估时就卡在这：后验分布没有解析形式，想要可信区间就必须从后验里采样，Matlab 自带的 randn、rand 那套常规武器根本没辙。搜索了一圈，所有答案都指向一个名字：马尔科夫链蒙特卡洛，也就是 MCMC。
MCMC 的基本思路相当反直觉：它不直接生成独立样本，而是构造一条马尔科夫链，让这条链的平稳分布恰好等于我们要采样的目标分布，然后让链自己跑足够长的时间，把链访问过的位置收集起来当样本。这个方法听起来像在绕远路，但它恰恰是贝叶斯推断、物理统计、计量金融等领域生成复杂分布数据的绝对主力。
这篇文章我打算用做项目的视角来讲，分三块内容：先解释 MCMC 为什么能干“只给密度函数就能采样”这件事，以及 Metropolis-Hastings 和 Gibbs 采样两种主流算法的工作原理；然后给出完整的 Matlab 实现，代码直接可以拿去跑，每一步的参数怎么设、为什么这么设都说明白；最后聊收敛诊断、样本质量评估和调参避坑——这些才是决定 MCMC 到底好不好用的关键。
不管你是正在做贝叶斯分析、需要合成数据做仿真实验，还是单纯想搞懂 MCMC 这两个词到底怎么回事，这篇文章都值得往下读。
1. 为什么“生成指定分布数据”比想象中难：MCMC 到底解决了什么问题
1.1 标准随机数生成器的天花板在哪里
如果目标分布恰好是均匀分布、正态分布、指数分布这类教科书分布，Matlab 确实一个命令就能搞定，rand、randn、exprnd 都现成。但真实项目的需求往往是另一回事：分布可能是多峰的、截断的、偏态的，甚至是只知道函数表达式、连积分都积不出来的怪东西。
常规替代方案我也都试过，各有各的死穴。逆变换采样要求目标分布的累积分布函数 CDF 能解析求逆，光这一条就把绝大多数复杂分布挡在门外。取舍采样（acceptance-rejection）倒是通用一些，但需要找到一根能完全罩住目标密度曲线的包络函数，在高维或多峰场景下，为了维持可接受的接受率，包络通常设计得极其宽松，采样效率低到没法工程化，尤其当维度一高，几乎等于不可用。
这些方法的天花板其实就一句话：它们大多要求目标分布的某种“全局信息”——CDF、包络、归一化常数——而真实问题常常只给得起“逐点函数值”。
1.2 只知道“未归一化的密度核”也能采样
这一点才是 MCMC 真正高明的地方。拿贝叶斯后验来举例：后验分布 π(θ|x) 正比于 似然 × 先验，而完整表达式的分母，也就是归一化常数，通常是一个没有解析解的高维积分。你算不出来这个积分，就画不出 CDF，就做不了逆变换采样。
但注意，MCMC 的接受准则里只需要两个位置的密度比值，归一化常数在分子分母里直接约掉了。这是什么概念？意味着面对一个密度函数，哪怕它的积分值我完全不知道，只要对任意一个输入 x 能算出一个正的、和真实密度成比例的数值，MCMC 就能照常采样。
我用一个琐碎的例子理解这件事：假设某个分布密度核长这样，p(x) ∝ exp(-x²/2) × (2 + sin(3x))，这玩意儿根本没有闭式积分，但你依然可以在任意点算出它的函数值。对 MCMC 来说这就够了。
1.3 马尔科夫链和蒙特卡洛为什么非要凑一块
“蒙特卡洛”好理解，就是用大量随机样本来逼近统计量，本质上是给积分做数值估计。可“马尔科夫链”扮演什么角色？
我习惯用一个逛商场的类比来解释。想象你逛一个商场，商场里不同区域人流热度不一样，有的地方好逛所以人多，有的地方冷清所以人少。你每走一步，转移倾向略微偏向热度高的区域：从冷清区往热闹区走的概率，比从热闹区往冷清区走的概率稍微大一点。只要这个“走下一步只看当前位置”的规则设计得合理，你逛足够久之后，停留在某个区域的时间比例，就会贴近那个区域的热度——也就是目标分布。
“下一步只和当前位置有关”就是马尔科夫的无记忆性质；“时间比例收敛到目标分布”就是链的平稳分布；“用这些逛出来的位置算均值、方差、分位数”就是蒙特卡洛那一半。马尔科夫链负责把链引到正确的分布上，蒙特卡洛负责把轨迹转换成统计结论。两者合体之后，数据生成问题就变成了“设计一条链”的问题。
2. 核心原理拆解：细致平衡条件与 Metropolis-Hastings 算法
2.1 细致平衡：链能稳定“逛”下去的数学保障
如果一条马尔科夫链有一个分布 π，使得从任意状态出发，跑一步之后新状态的边际分布还是 π，那就说 π 是这个链的平稳分布。真正让理论成立的是一个更强的充分条件——细致平衡：对任意两个状态 x 和 x'，都满足
π(x) P(x → x') = π(x') P(x' → x)
式子左边是“从 x 出发、一步走到 x'”的概率流量，右边是反向流量。细致平衡等于说，在平稳分布 π 之下，任意两状态之间的双向流量相等，净流量为零，所以链一旦进入这种状态就退不出来了。
构造 MCMC 算法的核心任务，就是把转移概率 P 设计成满足细致平衡。Metropolis 在 1953 年、Hastings 在 1970 年分别给出的方案，本质都是：先随便提议一个候选点，再用一个校正因子决定接不接受。
2.2 Metropolis-Hastings 的四步采样流程
标准的随机游走 Metropolis-Hastings（简称 MH）流程非常固定，一共四步：
从提议分布 q(x' | x_t) 生成一个候选点 x'。最常见的选择是给当前状态加一个高斯扰动：x' = x_t + σ·randn。
计算接受概率：
α = min(1, π(x')·q(x_t | x') / (π(x_t)·q(x' | x_t)))
如果提议分布对称，比如高斯分布，则 q(x_t | x') = q(x' | x_t)，接受概率直接简化成：
α = min(1, π(x') / π(x_t))
生成一个 [0,1] 均匀随机数 u，如果 u < α 就接受候选点，让链跳到 x'；否则停留在原状态 x_t。
把当前状态记入链中，回到第 1 步。
这里有个很容易困惑的地方：为什么目标密度更高的时候要接受，密度更低的时候偶尔也接受？答案是，如果链条只能向密度高的方向走，它迟早会困在某个局部峰顶永远出不来。允许以一定概率接受“下山”的移动，链才会具备探索整个分布空间的能力；而拒绝则保持原样，恰好让反向流量补齐，细致平衡才成立。所以接受率不是越高越好，这也是后面要重点调的地方。
2.3 一段可以直接跑的 Matlab 实现
先给一个最基础但完整的实现，目标分布选一个混合高斯：0.7 的 N(0,1) 叠加 0.3 的 N(5,2²)。代码用匿名函数定义目标密度，循环部分就是上面四步的直译。
% ===== Metropolis-Hastings 采样混合高斯分布 ===== rng(2026); % 固定随机种子，保证可复现 target_pdf = @(x) 0.7 * normpdf(x, 0, 1) + 0.3 * normpdf(x, 5, 2); N = 30000; % 总迭代次数 burn_in = 5000; % 燃烧期（前段样本要丢掉） sigma = 2.5; % 提议分布的标准差，也就是随机游走的步长 x_t = -3; % 初始状态 chain = zeros(N, 1); acc_num = 0; for i = 1:N x_star = x_t + sigma * randn(); % 1. 生成候选点 alpha = min(1, target_pdf(x_star) / target_pdf(x_t)); % 2. 接受概率 u = rand(); if u < alpha % 3. 是否接受 x_t = x_star; acc_num = acc_num + 1; end chain(i) = x_t; % 4. 记录当前状态 end samples = chain(burn_in+1:end); % 去掉燃烧期 accept_rate = acc_num / N; fprintf('接受率: %.2f%%\n', 100 * accept_rate); fprintf('样本均值: %.4f（理论值 1.5）\n', mean(samples));跑完之后样本均值应该接近 1.5，也就是 0.7×0 + 0.3×5 = 1.5 的混合期望；接受率大概落在 40% 上下。代码里每一步都标了注释，对应上面的流程，改 target_pdf 这个函数句柄，就能采样任意你能写出函数值的分布。
提示：
target_pdf返回的值必须恒为正。如果目标分布是在对数域里计算的，记得最后取 exp 还原，否则接受率会变成恒 0 或恒 1，采样结果完全失真。
3. 实战：用 Matlab 生成混合高斯数据的完整过程
3.1 为什么拿混合高斯当测试靶子
做数据生成的第一个原则就是：先别急着上真实目标分布，找一个理论性质完全已知的分布来验证算法。混合高斯是个绝佳的靶子，原因有三个。第一，它有多峰结构，峰与峰之间有低谷，最能考验链的“翻山越岭”能力；第二，它的真实均值、方差都能手算出来，方便和采样结果对比；第三，单变量场景下画轨迹图、直方图都直观，问题一眼就能暴露。
这种先验证后上生产的习惯，让我在后面的项目里少走了很多弯路。你不想在调试算法的同时还要连带着怀疑目标密度函数本身是否写错。
3.2 完整代码与参数化设计
把上一节的代码扩展一下，加上可视化和统计验证，形成一套可直接复用的脚本：
% ===== MCMC 数据生成完整示例：MH 采样混合高斯 ===== rng(2026); target_pdf = @(x) 0.7 * normpdf(x, 0, 1) + 0.3 * normpdf(x, 5, 2); N = 30000; burn_in = 5000; sigma = 2.5; x_t = -3; chain = zeros(N, 1); acc_num = 0; for i = 1:N x_star = x_t + sigma * randn(); alpha = min(1, target_pdf(x_star) / target_pdf(x_t)); if rand() < alpha x_t = x_star; acc_num = acc_num + 1; end chain(i) = x_t; end samples = chain(burn_in+1:end); sample_c = samples(1:2:end); % 稀疏化：每隔一个取一个，降低自相关 % 输出基本统计量 fprintf('接受率: %.2f%%\n', 100 * acc_num / N); fprintf('均值: %.4f（理论 1.5）\n', mean(sample_c)); fprintf('标准差: %.4f\n', std(sample_c)); % 轨迹图 + 直方图对比 figure('Color','w'); subplot(2,1,1); plot(1:N, chain, 'LineWidth', 0.3); hold on; xline(burn_in, '--r', 'LineWidth', 1.2); xlabel('迭代次数'); ylabel('采样值'); title('马尔科夫链轨迹图'); subplot(2,1,2); histogram(sample_c, 60, 'Normalization', 'pdf'); hold on; xx = linspace(min(sample_c), max(sample_c), 500); plot(xx, target_pdf(xx), 'r-', 'LineWidth', 2); xlabel('x'); ylabel('概率密度'); legend('MCMC 样本直方图', '真实密度', 'Location', 'northwest');单看样本直方图，两个峰的位置和真实密度曲线会重合得相当好，这就是生成成功最直观的证据。
3.3 步长这个参数，决定了链是“蠕动”还是“原地踏步”
MH 算法里唯一需要重点调试的参数就是提议分布的标准差 σ。步长的大小直接决定采样效率，我实测下来不同 σ 的表现差异非常明显：
- σ = 0.1：接受率会冲到 90% 以上，但每一步只挪一点点，链要花几十万步才可能从左侧峰逛到右侧峰，没逛到的地方就是采样盲区。轨迹图看起来像一条紧紧蜷在一起的细毛虫。
- σ = 2.5：接受率大约 40%，链在两个峰之间来回切换，轨迹图像毛线团但能铺满整个支撑域。这是最理想的区间。
- σ = 15：接受率掉到 10% 以下，链经常被“大跳跃”提议击中低密度区而拒绝，轨迹图出现大量水平线段，也就是长时间原地不动。
- σ = 50：接受率趋近于零，链几乎钉死在初始点附近，这组数据基本报废。
所以业界常年流传的经验是：把步长调到接受率落在 20% 到 50% 之间。这个经验不是严格的数学结论，而是大量实践验证的好用区间。调的时候可以写个小循环，扫描几个 σ 值，打印各自接受率，再挑一个中等的就行。
4. 高维场景的主力：Gibbs 采样原理与 Matlab 实现
4.1 Gibbs 采样的核心哲学：把高维拆成一维
MH 算法在高维空间会遇到一个问题：维度越高，接受率衰减越严重，因为随机游走提议在一个方向上移动时，其他维度的位置不变，整体移动到高密度邻域的概率指数级下降。Gibbs 采样则换了一条路——它根本不提议、不拒绝，而是依次从每个变量的条件分布中直接采样。
假设要采样的二元分布是 p(x1, x2)，Gibbs 的工作方式是一轮一轮地：
- 从条件分布 p(x1 | x2) 采样一个新 x1；
- 再从条件分布 p(x2 | x1) 采样一个新 x2；
- 得到一个新的二维样本点。
因为每一步用的都是真实的满条件分布，所以接受率恒为 1，不需要调任何提议参数。代价是，你必须能写出每个变量的条件分布，并且能从条件分布中采样。这听起来限制很大，但贝叶斯模型里大量共轭结构都能满足这个条件，所以它在分层模型、隐变量模型里应用极广。
4.2 二元高斯条件下的 Matlab 实现
二元联合正态是一个最经典的 Gibbs 演练场。设相关系数为 ρ，边际方差为 1，条件分布有解析形式：
p(x1 | x2) = N(ρ·x2, 1 − ρ²)
p(x2 | x1) = N(ρ·x1, 1 − ρ²)
写成 Matlab 代码非常短：
% ===== Gibbs 采样：二维正态分布 ===== rng(2023); rho = 0.8; N = 20000; burn_in = 2000; x1 = 0; x2 = 0; samples = zeros(N, 2); for i = 1:N x1 = rho * x2 + sqrt(1 - rho^2) * randn(); % 条件采样 x1 | x2 x2 = rho * x1 + sqrt(1 - rho^2) * randn(); % 条件采样 x2 | x1 samples(i, :) = [x1, x2]; end final = samples(burn_in+1:end, :); % 画出散点图，观察椭圆形状和相关性 figure('Color','w'); scatter(final(:,1), final(:,2), 8, [0.2 0.3 0.7], '.'); xlabel('x1'); ylabel('x2'); title('Gibbs 采样得到的二维正态样本'); fprintf('估计相关系数: %.3f（理论值 %.1f）\n', corr(final(:,1), final(:,2)), rho);散点图会呈现明显的沿对角线拉长的椭圆，相关系数估计值和 0.8 很接近。这个例子虽然简单，但它说明了 Gibbs 最关键的能力：高维联合分布可以通过一轮一轮交替采样一维条件分布来生成。
4.3 Metropolis-Hastings 和 Gibbs 怎么选：一张表说清楚
| 对比维度 | Metropolis-Hastings | Gibbs 采样 |
|---|---|---|
| 目标分布要求 | 只要任意点密度比值可算 | 每个变量的满条件分布必须可采样 |
| 参数调试 | 需要精心调提议分布的步长 | 没有提议分布，基本无参数 |
| 接受率 | 通常控制在 20% 到 50% | 恒为 1 |
| 高维表现 | 维度升高后接受率急剧恶化 | 条件分布质量好，维度相对无感 |
| 主要风险 | 多峰时容易困在单个峰里 | 变量强相关时链收敛很慢 |
一句话总结我的选型习惯：能用 Gibbs 就用 Gibbs，因为它安稳、不用调参；推导不出条件分布、或者条件分布没法采样的时候，才退回 MH 用提议分布硬闯。
5. 生成完不等于完事：收敛诊断与数据质量检验
5.1 燃烧期和稀疏化：样本不是从头到尾都能用
MCMC 样本里噪音最大的来源，就是链条最开头的那一段。我设定初始状态 x0 = -3 的时候，链并没有立刻进入目标分布，它需要一段时间“摸索”进高密度区。这段时间产生的前几千个样本如果直接混进结果，均值会被严重拉偏。这就是燃烧期，代码里用 burn_in 参数把它切掉。
另一个问题是链条相邻样本强相关。在 σ = 2.5 的例子里，相邻几步之间的状态高度相似，如果直接把全部样本当作独立样本用，算出来的方差会偏小、“有效信息量”被高估。处理办法是稀疏化：每隔 k 个样本取一个。我上面代码里的samples(1:2:end)就是这个用途。别把稀疏化当成强制步骤，如果后续分析本来就接受相关性，那全量保存也没问题；但如果要紧耦合统计推断，稀疏化能省掉很多解释上的麻烦。
5.2 三条必看的诊断图
我每次跑完 MCMC 都会画三张图，缺一不可：
- 轨迹图：直接画出链在每个迭代的取值。健康的链应该像毛线团一样充满整个取值区间，并且没有明显的“爬升”或“漂移”趋势；如果开头有很长一段单调爬坡，说明燃烧期没切够；如果出现长时间平台，说明链卡住了。
- 自相关图：看相隔 k 步的样本相关性。好的链自相关应该快速衰减到零附近；如果衰减很慢，说明链“蠕动”严重，需要加大步长或做稀疏化。Matlab 里
autocorr函数可以直接用，如果没装计量经济学工具箱，可以用xcorr归一化后自己画。 - 样本直方图叠加理论密度：这是最终的验收标准。把采样直方图和目标分布的真实密度画在同一张图上，直观比对峰的个数、位置、高度是否对得上。对不上的地方就是问题所在。
figure('Color','w'); subplot(1,3,1); plot(samples, 'LineWidth', 0.3); title('轨迹图'); xlabel('迭代'); ylabel('x'); subplot(1,3,2); autocorr(samples, 50); title('自相关图'); subplot(1,3,3); histogram(samples, 60, 'Normalization', 'pdf'); hold on; plot(xx, target_pdf(xx), 'r-', 'LineWidth', 1.5); title('密度对比');5.3 定量指标：Gelman-Rubin 诊断和有效样本量
光看图容易骗自己，所以我还会跑一个定量诊断。Gelman-Rubin 的思路是：并行跑 M 条链，而且每条链的初始状态故意选得差别很大，然后比较链间方差和链内方差。
具体做法不复杂：对每条链计算各自的样本均值和方差，得到组内方差 W；再计算所有链的总均值、以及链间均值差异带来的方差 B。如果链收敛了，链间方差应该远小于链内方差，两者的组合比例 R-hat 应该接近 1。经验阈值是 R-hat < 1.1 就认为收敛合格。
Matlab 里手动实现也不难，核心代码就是几行统计量计算：
% 伪代码示意：M 条链，每条长度 N chain_data = zeros(M, N); % ... 用不同初始值分别跑 M 次，填充 chain_data ... within = mean(var(chain_data, 0, 2)); % 组内方差均值 all_mean = mean(chain_data(:)); between = N / (M - 1) * sum((mean(chain_data, 2) - all_mean).^2); % 组间方差 Rhat = sqrt(((N - 1) / N * within + between / N) / within);判断标准很简单：R-hat 小于 1.1，放心用；大于 1.1，加大迭代次数或者检查步长，别硬着头皮拿未收敛的链生成数据。
5.4 我在实际项目里踩过的四个坑
这些坑没有一个是算法理论书里会专门警告的，全是我拿真实数据熬出来的：
- 只跑一条链就下结论。单链轨迹图可能看着正常，但它可能只拜访了其中一个峰。多峰分布尤其危险，我用混合高斯测试时就见过单链只采到左边峰的情况。现在我的习惯是至少三条链、不同初始值并行跑，用 R-hat 把关。
- 忘记切燃烧期。最典型的血案发生在初始状态取到极低密度区时，前几千个样本像一条从山脚爬向山顶的线，不切掉它，均值会被带偏一截。
- 把强相关的样本当独立样本做区间估计。MCMC 样本的有效量永远小于样本量，自相关性越强，差距越大。在计算置信区间、标准差之前，先看自相关图，必要时做稀疏化，否则你算出来的不确定性是偏乐观的。
- 在未归一化密度和对数密度之间切换时算错。很多实际项目的目标分布是在对数域计算的，比如 log-posterior。写代码时一旦忘记在算接受率前把对数差还原成比值，接受率就会直接告警：要么全接受要么全拒绝。排查这个问题只需要打印一轮接受率，看它是不是极端值即可。
6. 抛开代码看场景：这批数据生成代码能用在哪些地方
6.1 三个真实场景
第一个是贝叶斯后验采样，这是 MCMC 最经典的主场。后验分布没有解析表达式时，用 MH 或 Gibbs 从后验里生成样本，再直接算样本均值、分位数，得到参数估计和可信区间。我当初建模型时就是这么干的，一步到位解决了不确定性量化的问题。
第二个是合成数据做仿真实验。比如我写了一个新的参数辨识算法，但手里没有带标注的真实数据，那就先用 MCMC 按已知参数生成一批模拟数据，跑一遍算法，看它能不能把原始参数“挖”回来。这比拿真实验数据测试干净得多，因为真值是你自己放的，对错一目了然。
第三个是数据扩增和再采样。当训练集样本不平衡，或者你需要从某个先验知识中补充样本时，只要能把领域约束写成密度函数，MCMC 就能从里面往外掏数据。当然它比直接复制样本要慢，但在“必须服从指定分布”的硬约束下，这是为数不多靠谱的生成方式。
6.2 代码工程化的几条建议
我在把 MCMC 代码从实验脚本改造成可复用模块时，总结了几条实操建议：
- 把目标密度封装成独立的函数句柄，像上面代码里的
target_pdf。这样切换目标分布时只改一行，不用动采样循环。 - 固定随机种子。在脚本第一行写
rng(2026)，让每次运行结果完全一致。MCMC 自带随机性，不固定种子，昨天和今天的结果对不上，排查问题时你会非常痛苦。 - 多链并行。手动循环跑 M 条链时，可以用
parfor并行加速。链与链之间本来就相互独立，是天然适合并行的负载。 - 保存完整链而不是只存统计量。后面要做 R-hat、自相关分析时，缺原始链数据会很被动。内存允许的情况下，建议把整个 chain 矩阵存成
.mat文件归档。 - 先跑通现成函数再自己写。Matlab 统计工具箱里的
mhsample、slicesample可以直接用来快速验证思路；但自己的算法框架里，我还是推荐自己实现一遍采样循环，因为你必须理解每个参数的含义才能正确调优。
6.3 觉得 MH 高维太慢？下一个进阶方向是 HMC
到这里，基于 MCMC 的数据生成在 Matlab 里已经可以落地了。但如果你跑到高维模型，比如几十个参数的后验采样时发现效率暴跌，那不是实现的问题，而是随机游走的天然短板——MH 在维数升高以后，每一步的有效位移会急剧减小。
这时候的升级路径是 HMC，也就是哈密顿蒙特卡洛。它的思路是给链引入“动量”，利用目标分布的对数密度梯度信息来引导移动，像把一个小球放在势能场里，让惯性带着它穿过低概率区域而非一路乱逛。它对目标分布的平滑性有要求，但一旦用起来，维度诅咒会显著缓解，这也是现代概率编程库普遍采用它的原因。
建议的学习顺序是：先把这篇的 MH 和 Gibbs 跑熟，建立对马尔科夫链的直觉，再去看 HMC 和 NUTS。台阶是渐进的，直接上手 NUTS 容易变成调参机器，反而不利于理解原理。
最后聊一点我现在已经固化的使用习惯。每写一套 MCMC 数据生成代码，我都会强制自己做两件事：第一，先用一个已知解析答案的分布验证算法，确认样本均值、标准差、相关系数全部对得上，再跑去处理真实任务。这个过程很枯燥，但它能过滤掉绝大部分笔误和逻辑错误。第二，从第一天起就固定随机种子，把rng(...)写在脚本最前面，并在实验记录里标上这个种子和参数表。MCMC 生成本质上是随机过程，种子不固定，结果就没法复现，后续排查问题和写报告的时候，这种“笨功夫”比任何高级技巧都管用。
