企业官网建设流程全解析

简介：基于MCMC马尔科夫-蒙特卡洛抽样的Matlab仿真资源，面向本硕博及科研教学人群，适合用于马尔科夫链蒙特卡洛抽样算法的编程学习、教学演示与算法验证。资源共5个文件，以Matlab源程序为主（3个.m），另含1个avi操作录像和1个txt说明文档，压缩包整体仅753KB，轻量紧凑、即下即用。代码层面包含主运行脚本、被调用的子函数以及MCMC抽样算法实现，可直接在Matlab 2021a及以上版本运行；操作录像展示从路径设置、脚本执行到结果观察的完整过程，txt说明则补充了主程序与子函数的运行关系及注意事项，学习者可按视频快速复现抽样过程。目前已积累2381人学习下载，适合希望通过源码加演示系统理解马尔科夫链、建议分布与接受-拒绝机制的初学者，也便于教师在课堂上作为实验案例使用。

1. 你以为MCMC是黑箱，其实它就是一条能收敛的马尔科夫链

当后验分布的归一化常数求不出来、解析解不存在时，贝叶斯推断就卡住了。MCMC把问题从“求积分”换成“抽样”：构造一条以目标分布为平稳分布的马尔科夫链，沿着链子走下去，轨迹的样本分布自然逼近那个求不出来的后验，期望值用样本均值替代。这条思路看着绕，工程上却非常可靠。matlab仿真里最常用的两个入口是Metropolis-Hastings（MH）和Gibbs采样器：前者只需要未归一化的密度表达，后者把联合后验拆成一组好采的条件分布。这篇内容面向需要用MATLAB把MCMC跑通、又不想只按工具箱黑箱调参的人，从平稳分布原理讲到接受率、提议方差、收敛诊断和贝叶斯回归实战，代码可直接复制运行。

2. MCMC的马尔科夫链机制：平稳分布、细致平衡与混合效率

2.1 从“求不出积分”到“抽得出样本”：蒙特卡洛在数据增强里的角色

后验推断的通用形式是算后验期望 E[f(θ)|y] = ∫f(θ)π(θ|y)dθ，其中π(θ|y) ∝ L(y|θ)π(θ)。分子好写，分母是那个归一化常数∫L(y|θ)π(θ)dθ，维度一高就积不出来。蒙特卡洛的想法是：从π(θ|y)直接抽取大量样本，均值收敛到上述期望。问题在于“直接抽”在很多分布上做不到——逆变换要求可积的CDF，拒绝采样在高维空间接受率指数级下降，重要性采样在重尾时方差爆炸。

MCMC换了个思路：构造一条转移链，让它在长期运行中“生成”目标分布的样本，不需要求归一化常数，也不需要直接采样。这条链的关键不在样本是否独立，而在转移核P(x→x′)能使目标分布π在每一步之后都保持不变，也就是不动点条件π = πP。工程上更常用的是一个充分条件叫detailed balance（细致平衡）：π(x)P(x→x′) = π(x′)P(x′→x)，满足它，π一定是该链的平稳分布。MH算法的全部工作，就是设计转移核让这个条件成立：给定当前x，先从建议分布q(x′|x)抽候选，再以接受率决定是否移动。

2.2 转移核与平稳分布：MH如何满足detailed balance

MH的接受率形式是 α(x, x′) = min(1, π(x′)q(x|x′) / (π(x)q(x′|x)))。如果提议分布对称，比如高斯随机游走提议x′ = x + N(0, σ²)，q(x|x′) = q(x′|x)可以约掉，接受率退化为min(1, π(x′)/π(x))，这就是Metropolis采样。这个简化是matlab仿真里最常用的结构：写目标分布时只需要定义联合密度而不用管归一化常数，因为比值里常数消掉了。这正是MCMC最大的工程红利——很多应用问题里写得出似然×先验，但归一化常数根本没闭式。

接受率高低成了转移核设计里的核心变量。接受率过高意味着新样本离当前点很近，链在局部徘徊，混合慢；过低意味着每步都在拒绝，链原地踏步。高维情况下理论最优接受率在0.234附近（Roberts等人对高斯目标的渐近结果），一维目标大约0.44。仿真中不必精确逼近这个数，只需监控接受率落在0.2~0.5区间，偏离太多就调整提议方差。混合效率直接决定后面看到的自相关强度。

2.2.1 接受率公式的对称性简化与数值风险

细致平衡是局部条件，它确保任意一步转移中从x到x′的期望流等于反向流。把两侧对x积分，会发现π在转移作用下不变。这个证明在MATLAB中不必显式写出，但它决定了接受率的结构：分子分母必须包含前向概率π(x′)q(x′|x)和反向概率π(x)q(x|x′)，二者比值失衡时通过随机接受修正。很多仿真发散的问题，追根溯源是写接受率时把q项当成常量约掉，破坏了详细平衡条件。

2.3 burn-in、混合与自相关：样本不是独立同分布的

从任意初始值出发，链需要一段“预热”才能进入平稳区域，这一段称为burn-in，样本要丢弃。进入平稳后，相邻样本高度相关——这是马尔科夫链的固有属性，不是bug。相关越强，链的有效样本量越小。量化上使用自相关长度τ，有效样本量ESS = N/τ，τ越大信息越少。仿真时用trace plot（横轴迭代次数，纵轴样本值）直接看混合状态：理想情况下链在两个模式之间频繁穿梭，而不是长时间困在一侧。三张图画完，MCMC仿真才算真的“可见”了。

3. 用MATLAB从零写Metropolis-Hastings采样器：对准双峰混合高斯

3.1 目标分布与提议分布的选型：为什么用随机游走高斯提议

双峰混合高斯是MCMC最常用的试验场：分布有两个峰值，如果链只在一个峰附近徘徊，trace plot一眼就能看出来。设定目标分布为π(x) = 0.6·N(x; -2, 0.8²) + 0.4·N(x; 3, 1.2²)。这个分布没有解析采样方法，但密度表达式简单。提议分布采用对称高斯随机游走：x′ = x + σ·randn()。对称性的好处是接受率里没有q项，实现简单；缺点是高维时固定方差难以匹配不同维度尺度，后面自适应MCMC章节会解决。在确定性仿真场景下，σ的取值就是最重要的调参旋钮。

先定义对数密度函数。这里有一个资深用户常踩的坑：不要直接用exp相加再取log，两个高斯分量强度差异大时，exp(w)会下溢为0，log变成-Inf。应该在log域计算，用log-sum-exp技巧：

function logp = log_target(x) % 双峰高斯混合的对数密度，返回标量 % 每个分量的log权重 + log高斯核，归一化常数在MH比值中抵消 logw = zeros(size(x, 1), 2); % 第一个分量：均值-2，标准差0.8，权重0.6 logw(:, 1) = log(0.6) - 0.5 * ((x + 2) / 0.8).^2; % 第二个分量：均值3，标准差1.2，权重0.4 logw(:, 2) = log(0.4) - 0.5 * ((x - 3) / 1.2).^2; % log-sum-exp：先减最大值再取exp，防止指数下溢 m = max(logw, [], 2); logp = m + log(sum(exp(logw - m), 2)); end

这段函数的关键在于log-sum-exp：把公共最大值m提出来，exp(logw - m)的值域被控制在(0,1]，log(sum)不会爆炸。均值、标准差、权重是目标分布的物理意义参数；而- log(0.8·sqrt(2π))这类常数可以放心丢掉，因为MH接受率是比值。

3.2 主循环代码与接受率：log域计算避免下溢

MH主循环如下。以x0=0为初始状态，前5000次作为burn-in丢弃，后20000次保留：

% 主循环参数 N = 20000; % 采样总数 burnin = 5000; % 预热步数 sigma = 1.5; % 提议标准差，随机游走步长 x = 0; % 初始状态 samples = zeros(N, 1); acc = 0; % 接受计数器 for t = 1:(N + burnin) % 高斯随机游走提议：对称分布，q项在MH中抵消 x_star = x + sigma * randn(); % log域接受率：接受概率 = exp(logp(x_star) - logp(x)) log_alpha = log_target(x_star) - log_target(x); if log(rand()) < log_alpha x = x_star; acc = acc + 1; end % 丢弃burn-in阶段，其余样本保存 if t > burnin samples(t - burnin) = x; end end acceptance_rate = acc / (N + burnin); fprintf('接受率: %.3f\n', acceptance_rate);

主线是“提议-接受/拒绝-记录”。log(rand()) < log_alpha把接受概率转化为一次均匀随机数与阈值的比较，避免了对alpha再取exp的数值风险。计数器放在if内部，分母包含burnin，因为预热阶段的接受率同样反映提议方差是否合适。运行结束后接受率约在0.3附近说明步长合理；低于0.1要减小sigma，高于0.6则增大sigma。

3.2.1 提议方差σ怎么设：三个档位与接受率对照

sigma的选择直接决定采样质量，经验上分成三个档位：

提示：自相关强不是采样失败，而是有效样本量不足的信号。此时优先调sigma，而不是盲目增加总迭代数。

3.3 采样结果的三个验证图：trace、直方图、自相关

画验证图的代码集中在一次figure里：

subplot(2, 2, 1); plot(1:N, samples); title('Trace Plot'); xlabel('迭代次数'); ylabel('x'); subplot(2, 2, 2); histogram(samples, 60, 'Normalization', 'pdf'); hold on; % 把真实混合分布画上去做对比 xs = linspace(-6, 7, 400); true_pdf = 0.6 * normpdf(xs, -2, 0.8) + 0.4 * normpdf(xs, 3, 1.2); plot(xs, true_pdf, 'r-', 'LineWidth', 2); legend('MCMC样本', '真实密度'); subplot(2, 2, 3); autocorr(samples, 50); % 自相关图，横轴lag，纵轴自相关 fprintf('样本均值: %.3f\n', mean(samples)); fprintf('样本方差: %.3f\n', var(samples));

3.3.1 三个图分别回答什么问题

trace plot回答“链收敛了吗”。收敛的链，纵向抖动幅度稳定，两个峰之间周期性跨越；若长时间只在一个峰附近，要么是burn-in不够，要么是sigma太小。直方图回答“样本分布对吗”。MCMC样本的直方图与真实密度红色曲线贴合，说明抽样系统正确；贴合但形状粗糙，说明样本数不够或自相关太高。自相关图回答“信息量真的够吗”。自相关从lag=0的1.0衰减到0的速度越快越好；若滞后30仍然高于0.5，有效样本量会显著低于名义样本量。

操作视频的演示里通常把这三张图并列放在同一个figure窗口，运行结束一眼确认收敛状态。下一步顺着变异方向走：当变量维度增多，靠手调sigma不现实，就要引入Gibbs采样和自动协方差调节。

4. 从Gibbs到自适应MCMC：样本效率、R-hat与Geweke诊断

4.1 Gibbs采样：用满条件分布更新每一个分量

当目标分布是多维联合分布，而每个分量的满条件分布有闭型时，Gibbs采样比MH更高效——不需要调节提议方差，也没有拒绝步。做法是把θ拆成(θ1, …, θd)，每一轮依次从p(θi|θ_{-i}, y)采样，替换当前值。这个循环保持每个条件分布为平稳分布，组合起来就采样了联合后验。工程上最常见的场景是共轭先验下的贝叶斯模型：正态-正态、正态-逆伽马组合都有闭式条件分布。

MATLAB中实现Gibbs的骨架如下：

% Gibbs采样骨架：三变量目标，条件分布假设均有闭型 theta = zeros(3, 1); n_iter = 10000; samples_gibbs = zeros(n_iter, 3); % 伪代码示意，实际条件分布需按模型代入 for iter = 1:n_iter % 依次更新每个分量，直接以更新的值作为后续分量条件 theta(1) = conditional_sample_1(theta(2), theta(3)); theta(2) = conditional_sample_2(theta(1), theta(3)); theta(3) = conditional_sample_3(theta(1), theta(2)); samples_gibbs(iter, :) = theta'; end

更新顺序有讲究：第1个分量总是使用最新的θ2、θ3，第2个分量使用刚更新的θ1和旧的θ3。这种“始终用新值”的扫描方式称为systematic scan，是MATLAB中条件更新矩阵运算的标准写法，能加快混合速度。

4.2 自适应提议：2.38²/d缩放与协方差学习期

Gibbs要求条件分布可采样，但很多后验并不满足。回到MH框架，在多维情况下把单变量随机游走升级为多元高斯提议：x′ = x + N(0, σ²·Σ)，其中Σ是目标分布的协方差估计，σ是缩放因子。常用的经验法则是σ = 2.38 / sqrt(d)，d为目标维度，这个缩放因子在高斯目标下渐近最优。MATLAB里通过chol分解生成多元高斯样本：

% 自适应提议的协方差学习过程 d = 2; Sigma_adapt = eye(d) * 0.02; % 初始小协方差 L_adapt = chol(Sigma_adapt, 'lower'); N_total = 8000; learn_period = 3000; % 学习期，之后提议固定 x = zeros(d, 1); samples_adapt = zeros(N_total, d); for t = 1:N_total % 用当前协方差矩阵的Cholesky因子生成多元高斯增量 increment = L_adapt * randn(d, 1); x_star = x + 2.38 / sqrt(d) * increment; log_alpha = log_target_vec(x_star) - log_target_vec(x); if log(rand()) < log_alpha x = x_star; end samples_adapt(t, :) = x'; % 仅在学习期内更新协方差估计，学习期结束后固定提议 if t < learn_period && mod(t, 50) == 0 Sigma_adapt = cov(samples_adapt(max(1, t-500):t, :)); Sigma_adapt = Sigma_adapt + 1e-8 * eye(d); % 正则化防奇异 L_adapt = chol(Sigma_adapt, 'lower'); end end

这里两个必调参数：学习期learn_period和正则项系数1e-8。学习期必须小于总迭代数，结束后提议协方差固定，否则链的非马尔科夫性破坏收敛理论。正则化防止协方差矩阵在样本少时奇异导致chol失败。注意cov()接收的是过去500个样本形成的窗口，滑动窗口比全历史更能捕捉局部协方差结构。

提示：自适应MCMC只是在学习期用过去的样本来“调参”，正式采样仍要求提议不变。仿真中若看到学习期之后接受率骤降，多半是learn_period设置过短，协方差还没收敛就锁定了。

4.3 收敛诊断三件套：R-hat、有效样本量与Geweke检验

把MCMC当黑箱用时收敛诊断是保险。R-hat需要至少两条独立链，比较链间方差与链内方差：

function Rhat = compute_rhat(chains) % chains: M条链 × N次迭代，矩阵 [M, N] = size(chains); % 链均值的方差 → 链间方差B chain_means = mean(chains, 2); overall_mean = mean(chain_means); B = N / (M - 1) * sum((chain_means - overall_mean).^2); % 每条链内方差的平均 → 链内方差W W = mean(var(chains, 0, 2)); % 总方差估计混合公式 var_est = (N - 1) / N * W + B / N; Rhat = sqrt(var_est / W); end

计算时注意var(x, 0, 2)的维度参数：0表示除以N-1，2表示按行计算。R-hat越接近1越好。有效样本量ESS针对单条链的自相关结构：

function ess = effective_sample_size(x) % 通过自相关估计自相关长度tau，再得到ESS n = numel(x); xc = xcorr(x - mean(x), 'coeff'); rho = xc(n:end); % 保留非负lag一侧 % 截断自相关尾部：噪声对长lag求和影响大 max_lag = min(1000, n - 1); tau = 1 + 2 * sum(rho(2:max_lag + 1)); ess = n / tau; end

自相关尾部噪声会让tau被高估，max_lag上限正是为了压制这一点。Geweke检验取前10%与后50%的样本均值比较，构造的z统计量落在[-1.96, 1.96]内表示未检测到收敛异常。三者的用法对照：

仿真发散时先看这三项，再决定是加长学习期、增大迭代数，还是回头检查接受率和提议方差。

5. 实战闭环：用MCMC拟合贝叶斯线性回归并输出后验可信区间

5.1 模型设定与条件后验推导

把前面的技术落到最常见的回归场景。模型为y = Xβ + ε，ε ∼ N(0, σ²I)。配共轭先验：β ∼ N(β0, Λ0⁻¹)，σ² ∼ IG(a0, b0)。该设定下条件后验有闭型：β|σ², y ∼ N( (XᵀX/σ² + Λ0)⁻¹(Xᵀy/σ² + Λ0β0), (XᵀX/σ² + Λ0)⁻¹ )；σ²|β, y ∼ IG( a0 + n/2, b0 + (y - Xβ)ᵀ(y - Xβ)/2 )。第二个条件分布只依赖残差，这是Gibbs能闭环的关键。

生成仿真数据：

rng(42); n = 200; X = [ones(n, 1), randn(n, 1), randn(n, 1)]; % 截距 + 两个自变量 beta_true = [0.5; -1.2; 0.8]; y = X * beta_true + 0.5 * randn(n, 1);

5.2 Gibbs循环的MATLAB实现：按条件更新β和σ²

% 共轭先验参数 d = size(X, 2); beta0 = zeros(d, 1); Lambda0 = eye(d) * 0.01; a0 = 0.01; b0 = 0.01; % 初始值 beta = beta0; sigma2 = 1.0; n_iter = 12000; burnin = 2000; beta_samples = zeros(n_iter - burnin, d); sigma2_samples = zeros(n_iter - burnin, 1); for iter = 1:n_iter % 1. 更新beta：多元正态，协方差矩阵来自X'X/σ² + Λ的逆 Xbar = X' * X / sigma2 + Lambda0; beta_mean = Xbar \ (X' * y / sigma2 + Lambda0 * beta0); Sigma_beta = inv(Xbar); Lb = chol(Sigma_beta, 'lower'); beta = beta_mean + Lb * randn(d, 1); % 2. 更新sigma²：逆伽马，注意gamrnd第二个参数是scale不是rate resid = y - X * beta; a_post = a0 + n / 2; b_post = b0 + 0.5 * (resid' * resid); sigma2 = 1 / gamrnd(a_post, 1 / b_post); if iter > burnin beta_samples(iter - burnin, :) = beta'; sigma2_samples(iter - burnin) = sigma2; end end

两个常见的坑都在这里：gamrnd(a, b)第二个参数是scale，逆伽马的rate参数需要取倒数，所以写1/b_post；beta_mean用矩阵左除Xbar \ (...)而不是inv(Xbar)*...，数值稳定性更好，尤其当Xbar接近病态时。

5.3 输出后验结果并与OLS对比

采样结束后直接计算后验统计量：

beta_mean_mcmc = mean(beta_samples, 1)'; beta_ci = quantile(beta_samples, [0.025, 0.975], 1)'; beta_ols = (X' * X) \ (X' * y); fprintf('OLS: %.3f %.3f %.3f\n', beta_ols); fprintf('MCMC: %.3f %.3f %.3f\n', beta_mean_mcmc); fprintf('95%% CI:\n'); disp(beta_ci);

输出对照关系大致如下表（实际运行会有小数位波动，真实系数应落在CI内）：

这里的核心收益是“输出完整后验而不只是点估计”：OLS给出点和标准误，MCMC给出分布，可以随时改置信水平、计算P(β1 < 0)这类问题。收尾时跑一下effective_sample_size(beta_samples(:,2))，把ESS卡在400以上再汇报结果；若低于阈值就增加迭代数，而不是缩短burn-in。整套流程对照操作视频逐帧确认：burn-in段恰好对应trace plot从初始值摆动到平稳区域的时刻，协方差学习期结束对应接受率由震荡转平缓的转折点，这两个视觉标志对齐了，你手上这套基于MCMC马尔科夫-蒙特卡洛抽样matlab仿真才算真正跑通。

本文还有配套的精品资源，点击获取