简介基于MCMC马尔科夫-蒙特卡洛抽样的matlab仿真资源面向本硕博学生及科研教学人员定位于解决抽样算法编程落地与学习验证问题。资源共5个文件、整体仅753KB包含3个m源文件、1个txt说明和1个avi操作录像其中Runme_MCMC.m为主程序MCMC_metropolis_single.m实现Metropolis抽样func_E_Ms_int_single.m为辅助函数操作录像则完整演示运行流程。已有2381人学习下载适合需要在matlab中快速上手马尔科夫链蒙特卡洛抽样的读者。配合录像实际运行Runme_M_.m可深入理解状态转移、接受拒绝机制及积分近似等关键点同时避开直接运行子函数导致的报错高效完成算法复现与实验拓展。1. 当数值积分在高维空间失效MCMC是最后一道防线去年我在做一个贝叶斯参数估计的活后验密度只能算到未归一化形式维度是12。用Matlab自带的integral2试到三维就卡得不行网格求积更是直接内存爆炸。后来换成MCMC抽样几分钟就出了一组可靠的期望估计。MCMC的思路很反直觉与其费力去逼近一个高维积分不如构造一条马尔科夫链让链的稳态分布恰好等于目标分布π(x)然后沿着链取样本用样本均值替代积分。这个压缩包里的MCMC_metropolis_single.m就是干这件事的Runme_MCMC.m则是把采样、积分、画图串起来的主脚本。下面我从Metropolis-Hastings的实现原理讲起再逐步拆到积分计算、运行排错和进阶提速确保你能照着把仿真跑起来。2. Metropolis-Hastings的Matlab实现提议、接受率与对称性简化2.1 马尔科夫链的平稳分布从哪里来构造MCMC的核心是让链的转移核P(x→x)满足细致平衡条件π(x)P(x→x) π(x)P(x→x)。一旦满足这个等式π就是这条链的平稳分布。Metropolis-Hastings的思路是任意给一个提议分布q(x|x)然后以接受率α接受候选点使得实际转移概率P(x→x) q(x|x)·α(x→x)代入细致平衡条件就能解出α的形式。接受率公式是α min(1, [π(x)q(x|x)] / [π(x)q(x|x)])。当q是对称分布即q(x|x)q(x|x)时公式退化为α min(1, π(x)/π(x))。这就是经典的Metropolis算法也是随机游走Metropolis最常见的形态。高斯提议分布天然对称所以你在很多Matlab教材里看到的MCMC例子都默认用高斯随机游走本资源中的MCMC_metropolis_single.m也很可能沿用了这一简化。2.2 一个可运行的MCMC_metropolis_single.m版本虽然压缩包里已经有一个MCMC_metropolis_single.m但为了让你理解其内部逻辑我写一个同函数签名的示例版本。常见做法是把目标分布的对数密度传入避免计算原始密度时的数值下溢或溢出。function [samples, accept_rate] MCMC_metropolis_single(x0, sigma, N, log_target) % x0: 初始点行向量 % sigma: 高斯提议分布的标准差标量或与x同长度的向量 % N: 要采集的样本数量 % log_target: 函数句柄输入行向量x返回对数密度值 % 输出samples: N x d 的样本矩阵accept_rate: 实际接受率 d numel(x0); samples zeros(N, d); x x0(:); % 确保是行向量 acc 0; for t 1:N % 从对称高斯提议分布采样 x_prop x sigma .* randn(1, d); % 计算接受率的对数形式 log_alpha log_target(x_prop) - log_target(x); if log_alpha 0 alpha 1; else alpha exp(log_alpha); end % 按alpha概率接受候选点 if rand alpha x x_prop; acc acc 1; end samples(t, :) x; end accept_rate acc / N; end这段代码的核心是用对数密度做差而不是直接计算π(x)/π(x)。因为很多目标分布是未归一化的指数族比如高斯后验的对数形式是二次函数直接算比值会丢精度对数差则非常稳定。参数sigma的取值直接影响接受率sigma太小链走得慢自相关高sigma太大候选点经常跑到低密度区接受率低。一般经验是把接受率控制在20%40%之间。2.3 提议分布的选择与接受率调试不是所有场景都适合高斯随机游走下面的表格列出几种常见提议分布及适用情况仿真包里用的是第一行的高斯随机游走但对其他变体也该心里有数。提议分布转移形式优点缺点典型场景高斯随机游走x x σ·N(0,I)实现简单、对称需要调σ高维易低效连续参数中等维度对称拉普拉斯x x τ·Laplace(0,1)重尾能远跳有小概率大步长破坏链目标分布重尾独立提议x ~ q(x)不依赖当前状态q难以匹配π时接受率低目标分布形式已知随机游走混合多种σ按概率切换兼顾局部与全局参数增加多尺度目标分布调试sigma时我一般先跑一小段链比如5000次打印accept_rate。低于0.15就把sigma缩小高于0.5就加大。也可以用自适应算法在采样过程中动态调整那部分我放到最后一章讲。注意在Matlab里用randn时sigma如果是标量就直接乘在随机数上如果是向量就用点乘两种方式代表各维度使用相同的步长还是不同的步长。3. 用MCMC算积分期望E_Ms_int_single.m 的参数逻辑与收敛性3.1 蒙特卡洛积分的基本思想MCMC采样的最终目的是计算某个函数关于目标分布的期望E[f(x)] ∫f(x)π(x)dx。π如果是后验分布f可以是参数本身、平方值或者任何你关心的物理量。蒙特卡洛的做法是用链上样本的算术平均去逼近期望E[f(x)] ≈ (1/N) Σ_{i1}^{N} f(x_i)但问题在于马尔科夫链的前几个样本可能还停留在初始点附近没进入稳态这部分样本如果混入计算会产生较大偏差。另外相邻样本之间存在自相关简单平均的方差会被低估。所以在真正用样本做积分之前必须处理两个问题丢弃burn-in和考虑自相关性。3.2 E_Ms_int_single.m 的实现思路从文件名看E_Ms_int_single.m 很可能就是“期望值计算单链”的函数。它接收MCMC生成的样本矩阵和一个被积函数句柄返回期望值。我写一个简化但等价于常见实现逻辑的函数function E E_Ms_int_single(samples, burnin, f_handle) % samples: MCMC生成的样本矩阵每行是一个d维样本 % burnin: 需要丢弃的预热样本数标量 % f_handle: 被积函数句柄输入一行样本输出一个标量 % 输出E: 期望值估计 N size(samples, 1); if burnin N error(burnin必须小于样本总数); end x_use samples(burnin 1 : end, :); M size(x_use, 1); vals zeros(M, 1); for i 1:M vals(i) f_handle(x_use(i, :)); end E mean(vals); end这段代码先检查burnin是否越界然后截取稳态样本逐行调用f_handle计算函数值最后用mean求平均。这里f_handle传入的是形如(x) x(1)^2 x(2)的匿名函数。注意MCMC样本是逐行存储的所以f_handle必须接受一行向量而不是列向量否则会出现维度不匹配的错误。另外如果样本量很大逐行循环在Matlab中不是最优解可以用arrayfun或者sum向量化但为了可读性这里保留循环。3.3 burn-in与间隔参数的选取参数典型取值影响如何判断总样本数N1e4~1e6方差随N增大而减小计算有效样本量ESSburn-in1000~5000消除初始点影响看trace plot是否平稳间隔thin5~20降低自相关看自相关函数ACF在E_Ms_int_single这个函数里缩包版本可能只接收samples和被积函数burn-in放在外部处理。运行Runme_MCMC.m时你可以先画一下样本轨迹图观察链是否在大约几百步后进入一个相对稳定的区域。那个拐点之前的样本就该被丢弃。如果样本序列的自相关太大可以在取出样本时每间隔k个点取一个比如samples(burnin1:thin:end, :)。但要注意间隔过大等于减少了有效样本量不是万能药。4. 跑通Runme_MCMC.m路径、版本与常见报错排查4.1 运行环境和路径最容易挂的地方压缩包里有个操作录像0023.avi但我还是要强调常被忽略的那条Matlab左侧的当前文件夹窗口必须是工程所在路径。很多人直接双击Runme_MCMC.m然后Matlab内部出现了找不到函数的错误就是因为当前文件夹在别的目录下。Matlab解析脚本依赖时只从当前文件夹和搜索路径里找工程目录内的子文件不会被自动加入路径。版本方面摘要里写明用Matlab 2021a或更高版本测试。我实测用R2023b运行正常如果用2020b或更老版本可能会碰到arguments块或mustBePositive这类较新语法导致的解析错误。如果你不想升级遇到语法报错时可以先查一下是否是旧版不支持把相应语法改成传统nargin校验即可。4.2 Runme_MCMC.m 的典型流程主脚本的运行逻辑通常是设置目标分布参数 → 调用MCMC_metropolis_single采样 → 调用E_Ms_int_single计算期望 → 画轨迹图和直方图。我按照这个逻辑写一个主脚本骨架你可以对照缩包中的Runme_MCMC.m来看%% Runme_MCMC.m - 主脚本 clear; clc; close all; rng(42); % 固定随机种子复现结果 % 1. 设置目标分布这里用一个双峰高斯混合作例子 mu1 [-2, 2]; mu2 [3, -1]; log_target (x) log(0.5 * mvnpdf(x, mu1, 0.8*eye(2)) ... 0.5 * mvnpdf(x, mu2, 0.6*eye(2))); % 2. MCMC采样 x0 [0, 0]; sigma 0.5; N 20000; [samples, acc] MCMC_metropolis_single(x0, sigma, N, log_target); fprintf(接受率: %.2f%%\n, acc * 100); % 3. 计算期望如 E[x1^2 x2^2] burnin 5000; E E_Ms_int_single(samples, burnin, (x) x(1)^2 x(2)^2); fprintf(期望值: %.4f\n, E); % 4. 可视化 figure; subplot(2,1,1); plot(samples(:,1), LineWidth, 0.5); title(x1轨迹); subplot(2,1,2); hist3(samples, [30 30]); title(样本分布);这里的rng(42)固定随机种子保证每次运行结果一致。如果你要跑真实场景把log_target换成你自己的分布即可。注意MCMC_metropolis_single 和 E_Ms_int_single 都必须放在当前文件夹或已加入路径的文件夹中否则Matlab会报“未定义函数或变量”。4.3 收敛性诊断怎么判定链已经进入稳态光看期望值和接受率还不够我每次跑完MCMC都会做三项诊断。第一是轨迹图。把样本序列画出来如果链在某个均值附近来回震荡没有明显漂移就认为进入稳态。如果看到长段缓慢爬坡说明拒绝率太高或链没融合。第二是自相关图。计算当前样本与滞后k个样本的相关系数Matlab里用autocorr(samples(:,1))即可。自相关衰减越快的链越高效衰减慢说明步长太小需要加大sigma。第三是多链对比。用两三个差异很大的初始点比如某一维度正负10分别跑链观察它们最终是否混到同一分布区间。这个在Matlab里可以用parallel多个循环跑也可以在普通循环中串行跑。下面是简单的多链Gelman-Rubin诊断代码chains cell(3,1); x0s [-10, 10; 0, 0; 10, -10]; % 三个初始点 for c 1:3 chains{c} MCMC_metropolis_single(x0s(c,:), sigma, 10000, log_target); end % 计算每条链的后验均值看是否接近 means cellfun((c) mean(c(burnin1:end,1)), chains); disp(means);Gelman-Rubin的准确实现要计算链间和链内方差但如果你看到三个均值差距在0.5个标准差内基本可以认为收敛了。顺便提一句动力学蒙特卡洛KMC是物理上模拟原子扩散的另一种方法虽然都带“蒙特卡洛”名字但原理完全不同别用KMC的思路来诊断这里的MCMC链。4.4 常见报错与排错清单错误现象可能原因解决方法未定义函数或变量Runme_MCMC当前文件夹不是工程路径用cd到解压目录或点击左侧文件列表进入矩阵维度必须一致目标分布函数返回了列向量确保log_target输入输出都是行向量无法打开MCMC_metropolis_single.m文件名拼写不一致或不在路径中检查文件是否存在避免大小写问题运行录像打不开缺少解码器用VLC播放或改用PotPlayer结果与示例不符随机种子未固定在Runme前加rng(0)5. MCMC采样进阶自适应提议与多链冷热启动技巧5.1 自适应提议协方差固定sigma的随机游走在高维场景下常会陷入两种尴尬某些维度变化平缓需要大步长另一些维度陡峭需要小步长。一个经典的改进是Roberts和Rosenthal的自适应MCMC在运行中用历史样本的协方差来调节提议分布尺度。我一般这样实现% 在MCMC循环内部每过 batch 步就更新一次sigma batch 100; scale 2.38 / sqrt(d); % 理论最优标度 for t 1:N if t burnin mod(t, batch) 0 cov_user cov(samples(1:t-1, :)); sigma scale * chol(cov_user 0.1*eye(d), lower); end % 其余逻辑与前面一致... end注意自适应只能在链收敛后才能做否则会把初始阶段的偏差固化进协方差。另外加了自适应后理论上的平稳性质会变实际应用中要对丢失渐近性质有所取舍。我通常是在正式采样前跑一批预采样估计好协方差再用固定协方差跑正式链。5.2 冷热启动与并行多链当目标分布高度多峰时单个链容易被困在一个峰里。一个实用技巧是并行温度退火同时跑若干条链第i条链的目标分布为π(x)^{1/T_i}T_i大于1的链相当于把分布变得平坦更容易跨峰。然后周期性地在不同温度的链之间交换位置让高温链探索全局低温链细化局部。在Matlab中可以用parfor并行跑parfor t 1:4 T_list [1, 2, 4, 8]; chain_t MCMC_with_temp(x0, sigma, N, log_target, T_list(t)); chains{t} chain_t; end最后只保留T1这条链的结果。如果parfor在你自己机器上启动失败先检查并行池是否启动成功或者退化成普通for循环输出完全一致但速度慢一些。温度序列的选取还有一个经验最高温度对应的分布近似均匀分布这样交换才能有效跨越峰谷。你可以先跑一次高斯混合模型测试观察T8那链是否在峰之间高频跳跃如果是则温度序列设计合理。本文还有配套的精品资源点击获取 实战MCMC：Matlab实现Metropolis-Hastings解决高维贝叶斯积分
简介基于MCMC马尔科夫-蒙特卡洛抽样的matlab仿真资源面向本硕博学生及科研教学人员定位于解决抽样算法编程落地与学习验证问题。资源共5个文件、整体仅753KB包含3个m源文件、1个txt说明和1个avi操作录像其中Runme_MCMC.m为主程序MCMC_metropolis_single.m实现Metropolis抽样func_E_Ms_int_single.m为辅助函数操作录像则完整演示运行流程。已有2381人学习下载适合需要在matlab中快速上手马尔科夫链蒙特卡洛抽样的读者。配合录像实际运行Runme_M_.m可深入理解状态转移、接受拒绝机制及积分近似等关键点同时避开直接运行子函数导致的报错高效完成算法复现与实验拓展。1. 当数值积分在高维空间失效MCMC是最后一道防线去年我在做一个贝叶斯参数估计的活后验密度只能算到未归一化形式维度是12。用Matlab自带的integral2试到三维就卡得不行网格求积更是直接内存爆炸。后来换成MCMC抽样几分钟就出了一组可靠的期望估计。MCMC的思路很反直觉与其费力去逼近一个高维积分不如构造一条马尔科夫链让链的稳态分布恰好等于目标分布π(x)然后沿着链取样本用样本均值替代积分。这个压缩包里的MCMC_metropolis_single.m就是干这件事的Runme_MCMC.m则是把采样、积分、画图串起来的主脚本。下面我从Metropolis-Hastings的实现原理讲起再逐步拆到积分计算、运行排错和进阶提速确保你能照着把仿真跑起来。2. Metropolis-Hastings的Matlab实现提议、接受率与对称性简化2.1 马尔科夫链的平稳分布从哪里来构造MCMC的核心是让链的转移核P(x→x)满足细致平衡条件π(x)P(x→x) π(x)P(x→x)。一旦满足这个等式π就是这条链的平稳分布。Metropolis-Hastings的思路是任意给一个提议分布q(x|x)然后以接受率α接受候选点使得实际转移概率P(x→x) q(x|x)·α(x→x)代入细致平衡条件就能解出α的形式。接受率公式是α min(1, [π(x)q(x|x)] / [π(x)q(x|x)])。当q是对称分布即q(x|x)q(x|x)时公式退化为α min(1, π(x)/π(x))。这就是经典的Metropolis算法也是随机游走Metropolis最常见的形态。高斯提议分布天然对称所以你在很多Matlab教材里看到的MCMC例子都默认用高斯随机游走本资源中的MCMC_metropolis_single.m也很可能沿用了这一简化。2.2 一个可运行的MCMC_metropolis_single.m版本虽然压缩包里已经有一个MCMC_metropolis_single.m但为了让你理解其内部逻辑我写一个同函数签名的示例版本。常见做法是把目标分布的对数密度传入避免计算原始密度时的数值下溢或溢出。function [samples, accept_rate] MCMC_metropolis_single(x0, sigma, N, log_target) % x0: 初始点行向量 % sigma: 高斯提议分布的标准差标量或与x同长度的向量 % N: 要采集的样本数量 % log_target: 函数句柄输入行向量x返回对数密度值 % 输出samples: N x d 的样本矩阵accept_rate: 实际接受率 d numel(x0); samples zeros(N, d); x x0(:); % 确保是行向量 acc 0; for t 1:N % 从对称高斯提议分布采样 x_prop x sigma .* randn(1, d); % 计算接受率的对数形式 log_alpha log_target(x_prop) - log_target(x); if log_alpha 0 alpha 1; else alpha exp(log_alpha); end % 按alpha概率接受候选点 if rand alpha x x_prop; acc acc 1; end samples(t, :) x; end accept_rate acc / N; end这段代码的核心是用对数密度做差而不是直接计算π(x)/π(x)。因为很多目标分布是未归一化的指数族比如高斯后验的对数形式是二次函数直接算比值会丢精度对数差则非常稳定。参数sigma的取值直接影响接受率sigma太小链走得慢自相关高sigma太大候选点经常跑到低密度区接受率低。一般经验是把接受率控制在20%40%之间。2.3 提议分布的选择与接受率调试不是所有场景都适合高斯随机游走下面的表格列出几种常见提议分布及适用情况仿真包里用的是第一行的高斯随机游走但对其他变体也该心里有数。提议分布转移形式优点缺点典型场景高斯随机游走x x σ·N(0,I)实现简单、对称需要调σ高维易低效连续参数中等维度对称拉普拉斯x x τ·Laplace(0,1)重尾能远跳有小概率大步长破坏链目标分布重尾独立提议x ~ q(x)不依赖当前状态q难以匹配π时接受率低目标分布形式已知随机游走混合多种σ按概率切换兼顾局部与全局参数增加多尺度目标分布调试sigma时我一般先跑一小段链比如5000次打印accept_rate。低于0.15就把sigma缩小高于0.5就加大。也可以用自适应算法在采样过程中动态调整那部分我放到最后一章讲。注意在Matlab里用randn时sigma如果是标量就直接乘在随机数上如果是向量就用点乘两种方式代表各维度使用相同的步长还是不同的步长。3. 用MCMC算积分期望E_Ms_int_single.m 的参数逻辑与收敛性3.1 蒙特卡洛积分的基本思想MCMC采样的最终目的是计算某个函数关于目标分布的期望E[f(x)] ∫f(x)π(x)dx。π如果是后验分布f可以是参数本身、平方值或者任何你关心的物理量。蒙特卡洛的做法是用链上样本的算术平均去逼近期望E[f(x)] ≈ (1/N) Σ_{i1}^{N} f(x_i)但问题在于马尔科夫链的前几个样本可能还停留在初始点附近没进入稳态这部分样本如果混入计算会产生较大偏差。另外相邻样本之间存在自相关简单平均的方差会被低估。所以在真正用样本做积分之前必须处理两个问题丢弃burn-in和考虑自相关性。3.2 E_Ms_int_single.m 的实现思路从文件名看E_Ms_int_single.m 很可能就是“期望值计算单链”的函数。它接收MCMC生成的样本矩阵和一个被积函数句柄返回期望值。我写一个简化但等价于常见实现逻辑的函数function E E_Ms_int_single(samples, burnin, f_handle) % samples: MCMC生成的样本矩阵每行是一个d维样本 % burnin: 需要丢弃的预热样本数标量 % f_handle: 被积函数句柄输入一行样本输出一个标量 % 输出E: 期望值估计 N size(samples, 1); if burnin N error(burnin必须小于样本总数); end x_use samples(burnin 1 : end, :); M size(x_use, 1); vals zeros(M, 1); for i 1:M vals(i) f_handle(x_use(i, :)); end E mean(vals); end这段代码先检查burnin是否越界然后截取稳态样本逐行调用f_handle计算函数值最后用mean求平均。这里f_handle传入的是形如(x) x(1)^2 x(2)的匿名函数。注意MCMC样本是逐行存储的所以f_handle必须接受一行向量而不是列向量否则会出现维度不匹配的错误。另外如果样本量很大逐行循环在Matlab中不是最优解可以用arrayfun或者sum向量化但为了可读性这里保留循环。3.3 burn-in与间隔参数的选取参数典型取值影响如何判断总样本数N1e4~1e6方差随N增大而减小计算有效样本量ESSburn-in1000~5000消除初始点影响看trace plot是否平稳间隔thin5~20降低自相关看自相关函数ACF在E_Ms_int_single这个函数里缩包版本可能只接收samples和被积函数burn-in放在外部处理。运行Runme_MCMC.m时你可以先画一下样本轨迹图观察链是否在大约几百步后进入一个相对稳定的区域。那个拐点之前的样本就该被丢弃。如果样本序列的自相关太大可以在取出样本时每间隔k个点取一个比如samples(burnin1:thin:end, :)。但要注意间隔过大等于减少了有效样本量不是万能药。4. 跑通Runme_MCMC.m路径、版本与常见报错排查4.1 运行环境和路径最容易挂的地方压缩包里有个操作录像0023.avi但我还是要强调常被忽略的那条Matlab左侧的当前文件夹窗口必须是工程所在路径。很多人直接双击Runme_MCMC.m然后Matlab内部出现了找不到函数的错误就是因为当前文件夹在别的目录下。Matlab解析脚本依赖时只从当前文件夹和搜索路径里找工程目录内的子文件不会被自动加入路径。版本方面摘要里写明用Matlab 2021a或更高版本测试。我实测用R2023b运行正常如果用2020b或更老版本可能会碰到arguments块或mustBePositive这类较新语法导致的解析错误。如果你不想升级遇到语法报错时可以先查一下是否是旧版不支持把相应语法改成传统nargin校验即可。4.2 Runme_MCMC.m 的典型流程主脚本的运行逻辑通常是设置目标分布参数 → 调用MCMC_metropolis_single采样 → 调用E_Ms_int_single计算期望 → 画轨迹图和直方图。我按照这个逻辑写一个主脚本骨架你可以对照缩包中的Runme_MCMC.m来看%% Runme_MCMC.m - 主脚本 clear; clc; close all; rng(42); % 固定随机种子复现结果 % 1. 设置目标分布这里用一个双峰高斯混合作例子 mu1 [-2, 2]; mu2 [3, -1]; log_target (x) log(0.5 * mvnpdf(x, mu1, 0.8*eye(2)) ... 0.5 * mvnpdf(x, mu2, 0.6*eye(2))); % 2. MCMC采样 x0 [0, 0]; sigma 0.5; N 20000; [samples, acc] MCMC_metropolis_single(x0, sigma, N, log_target); fprintf(接受率: %.2f%%\n, acc * 100); % 3. 计算期望如 E[x1^2 x2^2] burnin 5000; E E_Ms_int_single(samples, burnin, (x) x(1)^2 x(2)^2); fprintf(期望值: %.4f\n, E); % 4. 可视化 figure; subplot(2,1,1); plot(samples(:,1), LineWidth, 0.5); title(x1轨迹); subplot(2,1,2); hist3(samples, [30 30]); title(样本分布);这里的rng(42)固定随机种子保证每次运行结果一致。如果你要跑真实场景把log_target换成你自己的分布即可。注意MCMC_metropolis_single 和 E_Ms_int_single 都必须放在当前文件夹或已加入路径的文件夹中否则Matlab会报“未定义函数或变量”。4.3 收敛性诊断怎么判定链已经进入稳态光看期望值和接受率还不够我每次跑完MCMC都会做三项诊断。第一是轨迹图。把样本序列画出来如果链在某个均值附近来回震荡没有明显漂移就认为进入稳态。如果看到长段缓慢爬坡说明拒绝率太高或链没融合。第二是自相关图。计算当前样本与滞后k个样本的相关系数Matlab里用autocorr(samples(:,1))即可。自相关衰减越快的链越高效衰减慢说明步长太小需要加大sigma。第三是多链对比。用两三个差异很大的初始点比如某一维度正负10分别跑链观察它们最终是否混到同一分布区间。这个在Matlab里可以用parallel多个循环跑也可以在普通循环中串行跑。下面是简单的多链Gelman-Rubin诊断代码chains cell(3,1); x0s [-10, 10; 0, 0; 10, -10]; % 三个初始点 for c 1:3 chains{c} MCMC_metropolis_single(x0s(c,:), sigma, 10000, log_target); end % 计算每条链的后验均值看是否接近 means cellfun((c) mean(c(burnin1:end,1)), chains); disp(means);Gelman-Rubin的准确实现要计算链间和链内方差但如果你看到三个均值差距在0.5个标准差内基本可以认为收敛了。顺便提一句动力学蒙特卡洛KMC是物理上模拟原子扩散的另一种方法虽然都带“蒙特卡洛”名字但原理完全不同别用KMC的思路来诊断这里的MCMC链。4.4 常见报错与排错清单错误现象可能原因解决方法未定义函数或变量Runme_MCMC当前文件夹不是工程路径用cd到解压目录或点击左侧文件列表进入矩阵维度必须一致目标分布函数返回了列向量确保log_target输入输出都是行向量无法打开MCMC_metropolis_single.m文件名拼写不一致或不在路径中检查文件是否存在避免大小写问题运行录像打不开缺少解码器用VLC播放或改用PotPlayer结果与示例不符随机种子未固定在Runme前加rng(0)5. MCMC采样进阶自适应提议与多链冷热启动技巧5.1 自适应提议协方差固定sigma的随机游走在高维场景下常会陷入两种尴尬某些维度变化平缓需要大步长另一些维度陡峭需要小步长。一个经典的改进是Roberts和Rosenthal的自适应MCMC在运行中用历史样本的协方差来调节提议分布尺度。我一般这样实现% 在MCMC循环内部每过 batch 步就更新一次sigma batch 100; scale 2.38 / sqrt(d); % 理论最优标度 for t 1:N if t burnin mod(t, batch) 0 cov_user cov(samples(1:t-1, :)); sigma scale * chol(cov_user 0.1*eye(d), lower); end % 其余逻辑与前面一致... end注意自适应只能在链收敛后才能做否则会把初始阶段的偏差固化进协方差。另外加了自适应后理论上的平稳性质会变实际应用中要对丢失渐近性质有所取舍。我通常是在正式采样前跑一批预采样估计好协方差再用固定协方差跑正式链。5.2 冷热启动与并行多链当目标分布高度多峰时单个链容易被困在一个峰里。一个实用技巧是并行温度退火同时跑若干条链第i条链的目标分布为π(x)^{1/T_i}T_i大于1的链相当于把分布变得平坦更容易跨峰。然后周期性地在不同温度的链之间交换位置让高温链探索全局低温链细化局部。在Matlab中可以用parfor并行跑parfor t 1:4 T_list [1, 2, 4, 8]; chain_t MCMC_with_temp(x0, sigma, N, log_target, T_list(t)); chains{t} chain_t; end最后只保留T1这条链的结果。如果parfor在你自己机器上启动失败先检查并行池是否启动成功或者退化成普通for循环输出完全一致但速度慢一些。温度序列的选取还有一个经验最高温度对应的分布近似均匀分布这样交换才能有效跨越峰谷。你可以先跑一次高斯混合模型测试观察T8那链是否在峰之间高频跳跃如果是则温度序列设计合理。本文还有配套的精品资源点击获取