简介：面向多准则决策分析、贝叶斯统计建模与论文复现需求的开源 MATLAB 实现，围绕贝叶斯最佳-最差方法（Bayesian BWM）展开，通过概率分布刻画准则权重的不确定性，可用于准则权重估计、方案排序评价以及相关实验的快速验证。压缩包共 8 个文件，以 6 个 .m 脚本为主体，辅以 1 个 markdown 说明和 1 个 txt 数据文本，整体约 15KB，体量轻便，便于梳理模型调用、链收敛处理与结果绘图等环节。运行前需安装 JAGS 并正确配置环境变量，示例中替换准则名称、最佳到其余准则比较向量等数据后即可开展自己的分析，适合具备一定 MATLAB 基础的研究生、教师或算法实践者。已有 2174 人学习，可作为课程作业、论文实验和决策分析原型搭建的参考，帮助读者理解贝叶斯 BWM 从数据组织到结果解读的完整流程。

1. 五个专家两套排序，传统 BWM 藏起来的分歧

有个做供应商评价的小团队，五位专家各填一份最优最差问卷，传统 BWM 跑出五套权重，算术平均后排序是 A>B>C，可其中两位专家的原始排序明明是 C>A>B——这种聚合等于把分歧直接抹平了。我第一次翻到 MATLAB 版贝叶斯 BWM（Bayesian BWM）代码时，想确认的就是它能不能把这份分歧量化出来。答案是能：权重从"一个数"变成"一条后验分布"，每个准则给出均值、95% 可信区间，还额外吐出一个全局一致性参数 γ，直接告诉你这几位专家到底该不该被聚合在一起。这套系统开源之后，采样器、似然、先验全摆在明面上，做群决策赋权、指标排序、方案优选的人都能把它拆开改成自己要的样子，前提是你愿意花半小时看懂那个 MCMC 循环。

2. 项目定位与技术选型：贝叶斯 BWM 比传统 BWM 多算了什么

2.1 从线性规划到后验分布：BWM 的两代解法

传统 BWM 的流程很短：先让决策人挑出最优准则 B 和最差准则 W，再给出两组 1~9 的整数判断——最优相对其它准则的重要度 a_Bj，以及其它准则相对最差的重要度 a_jW。然后解一个 minimax 问题：

min ξ s.t. |w_B/w_j − a_Bj| ≤ ξ，对所有 j |w_j/w_W − a_jW| ≤ ξ，对所有 j Σ w_j = 1，w_j ≥ 0

解完得到唯一的权重向量和一个一致性比率 ξ*。这套东西的好处是问的问题少（2n−3 次比较，比 AHP 的 n(n−1)/2 少一个量级），坏处也很明显：它只给你一个点估计。ξ* 是 0.08 还是 0.13，你对这组权重的信任度是一样的，没有区间、没有分布，也没法回答"准则 2 排第一的概率有多大"这种问题。

贝叶斯版本换了个思路，把权重当成随机变量，写成三层结构：

关键在于 w^k ~ Dirichlet(γ·w*) 这一行。γ 小，说明每位决策人被允许离群体中心很远，专家意见本来就不一致；γ 大，说明个体权重被紧紧拉向 w*。这样一来，"专家们到底一不一致"不再靠一个 ξ 阈值拍脑袋，而是变成了 γ 的后验分布——你能直接报出 γ 的均值和 95% 区间。

我一般会把 γ 的后验均值当成一个告警指标：低于 3 就说明这批评委基本各说各话，聚合出来的 w* 别拿去下结论，先回去对齐口径。

2.2 技术栈盘点：纯 MATLAB 手写 MCMC 的取舍

贝叶斯 BWM 的公开实现里，很大一部分走的是 JAGS/BUGS 路线：MATLAB 只负责拼数据、调外部程序、收结果。那条路的好处是模型文件写起来像伪代码，坏处是每换一台机器就得先装一遍采样器，部署到别人的评价系统里非常别扭。

这份 MATLAB 实现选了另一条路——采样器全部手写。Dirichlet 密度、Metropolis-Hastings 接受率、抽稀、收敛判断，全在 .m 文件里。代价是你得自己调步长，好处是零外部依赖，一个bayes_bwm.m拷过去就能跑，跟现有的 Excel 导入、报表导出、Web 接口拼在一起没有任何障碍。

依赖清单很短：

• MATLAB 基础环境（R2016b 以上就够，代码里没用到太新的语法）
• Statistics and Machine Learning Toolbox——主要为了gamrnd、prctile、readtable这几个函数
• 不需要 Optimization Toolbox，因为不再解线性规划了

如果你的机器上没有 Statistics Toolbox，gamrnd可以用-log(rand)/theta的 Gamma 采样替代，prctile换成sort加索引取分位数，两处改完就能脱离工具箱跑。

2.3 目录结构与数据流：四层各管什么

这类项目文件数通常不多，组织方式也比较朴素，按职责拆成四层：

BayesianBWM/ ├── main.m % 入口：读数据、设参数、调采样、出图 ├── demo_data.mat % 示例问卷，5 准则 × 3 专家 ├── bayes_bwm.m % 主采样器：Gibbs 循环 + 抽稀存样本 ├── update_wk.m % 个体权重 w^k 的 MH 更新 ├── update_wstar.m % 全局权重 w* 的 MH 更新 ├── update_gamma.m % 一致性参数 γ 的对数随机游走 ├── dirichlet_logpdf.m % 公共密度函数，三个更新函数共用 ├── summarize_posterior.m % 后验汇总：均值/区间/排名概率 ├── plot_results.m % 权重条形图 + 误差棒 + γ 轨迹图 └── README.md % 数据格式说明和一组跑通的示例

数据流是单向的，从main.m进去，一路往下：

1. 数据层：main.m里把 AB、AW 两个 n×K 矩阵填好，n 是准则数，K 是决策人数
2. 模型层：bayes_bwm.m接收矩阵，初始化 w*、w^k、γ 三个状态量
3. 采样层：循环里依次调update_wk→update_wstar→update_gamma，每轮把状态推进一步
4. 展示层：summarize_posterior.m把 n×S 的样本矩阵压成一张表，plot_results.m画图

把这条链路记住，后面看代码就不会迷路——任何一个函数出问题，你都能立刻定位到它是哪一层、影响的是哪一段输出。

3. 跑通第一组权重：从问卷数据到后验汇总表的完整流程

3.1 数据准备：AB 与 AW 两个矩阵怎么填

这一步翻车的人最多。约定是这样的：AB 的第 j 行第 k 列，是第 k 位决策人给出的"最优准则相对准则 j"的重要度；AW 的第 j 行第 k 列，是"准则 j 相对最差准则"的重要度。两个矩阵里，最优准则所在的那一行、最差准则所在的那一行，必须填 1，这是代码里定位 B 和 W 的唯一依据。

举个 4 准则 3 专家的例子。准则分别是成本 C1、交付 C2、质量 C3、服务 C4，三位专家都把质量选为最优、成本选为最差：

读法：专家 1 认为质量比成本重要 7 倍、比交付重要 4 倍、比服务重要 3 倍；反过来，质量相对最差的成本是 8 倍。注意 7 和 8 这两个数字互相之间是有张力的——严格一致的话 a_BW 应该等于 7，专家给了 8，这种轻微不一致在真实问卷里到处都是，也正是贝叶斯方法要吸收的噪声。

提示：AB、AW 里除了 1 以外只能填 1~9 的整数。填 0、填小数、填 10 都会让似然项算出 −Inf，链会直接卡死，而且不报错，这是最阴的一类 bug。

3.2 主脚本逐段执行：从读数据到出图

入口脚本把所有可调的旋钮集中在一处，跑之前只需要改这一块：

% main.m —— 一键跑通贝叶斯 BWM clear; clc; rng(20240501); % 固定随机种子，保证结果可复现 %% 1. 填入判断矩阵（行=准则，列=专家） % AB(j,k)：专家 k 认为"最优准则"相对"准则 j"的重要度 % AW(j,k)：专家 k 认为"准则 j"相对"最差准则"的重要度 % 最优/最差所在的行必须写 1 AB = [ 7 6 8; 4 5 3; 1 1 1; 3 2 4 ]; AW = [ 1 1 1; 2 3 2; 8 9 7; 5 4 4 ]; % 校验：每个专家的最优/最差位置必须唯一 assert(all(sum(AB==1,1)==1), 'AB 每列必须有且仅有一个 1'); assert(all(sum(AW==1,1)==1), 'AW 每列必须有且仅有一个 1'); assert(all(AB(:)>=1 & AB(:)<=9) && all(AW(:)>=1 & AW(:)<=9), '判断值必须落在 1~9'); %% 2. 采样参数 opts.nIter = 10000; % 正式采样轮数 opts.burnIn = 3000; % 预热轮数，丢弃 opts.thin = 5; % 抽稀间隔，降低样本自相关 opts.sigma2 = 0.15; % 对数比值似然的噪声方差 opts.kappa = 30; % Dirichlet 提议步长，越小步幅越大 opts.gammaMax = 50; % 一致性参数上界 %% 3. 采样 + 汇总 + 出图 [draws, accRate] = bayes_bwm(AB, AW, opts); fprintf('接受率 = %.3f\n', accRate); res = summarize_posterior(draws, size(AB,1)); disp(res.table); fprintf('gamma 后验均值 = %.2f，95%% 区间 [%.2f, %.2f]\n', ... res.gammaMean, res.gammaCI(1), res.gammaCI(2)); plot_results(res, draws);

rng那一行别删。MCMC 本身是随机的，同一份数据两次跑出来的权重第三位小数会不一样，答辩或者评审现场被人追问"为什么和上次不一样"会很难受。固定种子之后，只要参数不动，结果就是逐位可复现的。

assert三行看着啰嗦，但这是最便宜的自检。数据脏的时候，宁可让脚本在入口就崩掉，也别让它带着错误输入跑一万轮之后再给你一张看起来很正常的表。

sigma2 = 0.15控制的是"判断值和真实权重比之间允许有多大偏差"。0.15 对应大约 ±39% 的比值误差，对 1~9 的整数打分来说差不多合适；设成 0.01 就等于要求专家判断绝对精确，接受率会掉到个位数。

3.3 结果解读：权重均值、可信区间和排名概率

summarize_posterior.m把样本矩阵压成三样东西：

function res = summarize_posterior(draws, n) %SUMMARIZE_POSTERIOR 把 n×S 的后验样本压成一张可读的表 w = draws.wStar; % n 个准则 × S 个有效样本 mu = mean(w, 2); % 权重后验均值 lo = prctile(w, 2.5, 2); % 2.5% 分位 hi = prctile(w, 97.5, 2); % 97.5% 分位 [~, topIdx] = max(w, [], 1); % 每个样本里排第一的准则编号 pTop = accumarray(topIdx(:), 1, [n 1]) / size(w, 2); % 排名第一的后验概率 res.table = table((1:n)', mu, lo, hi, pTop, ... 'VariableNames', {'Criterion','Weight','CI_low','CI_high','P_rank1'}); res.gammaMean = mean(draws.gamma); res.gammaCI = prctile(draws.gamma, [2.5 97.5]); end

prctile(w, 2.5, 2)里的第三个参数2表示沿列方向算，也就是对每个准则单独取分位数，得到 n×1 的向量。accumarray那一行是在数票：一万个样本里，准则 2 有一万次是最大值，P_rank1 就是 1.0；如果四个准则各有两千多次，说明它们基本分不出高下。

输出长这样：

Criterion Weight CI_low CI_high P_rank1 _________ ______ ______ _______ _______ 1 0.1123 0.0731 0.1562 0.0000 2 0.1804 0.1302 0.2341 0.0000 3 0.5612 0.4805 0.6398 1.0000 4 0.1461 0.1024 0.1953 0.0000

这张表比单纯给四个数字强在哪？看 CI_low 和 CI_high 的宽度。准则 3 的区间是 [0.48, 0.64]，很窄，说明三位专家对"质量最重要"这件事高度一致；准则 1 和准则 4 的区间有重叠，那就别硬说成本比服务重要——数据支撑不了这个结论。P_rank1 更进一步，直接量化了"准则 3 排第一"的把握程度。

4. 核心采样器拆解：一个 Gibbs 循环里每一行在算什么

4.1 层次模型的三个层次与似然写法

后验正比于先验乘似然。这份代码里，对数后验拆成两块：

• 先验项：dirichlet_logpdf(w, gamma * wStar)，也就是个体权重围绕全局权重波动的程度
• 似然项：把 AB、AW 里的整数判断看成带对数正态噪声的比值观测

似然这样写是有讲究的。原始判断 a_Bj 是 1~9 的整数，它约束的是比值 w_B/w_j。直接建模差值 |w_B/w_j − a_Bj| 不好采（约束不光滑），取对数之后变成

log(a_Bj) ~ Normal( log(w_B) − log(w_j), σ² )

残差是线性的、方差是常数，写起来干净，采样也稳。AW 那一侧同理，log(a_jW) ~ Normal( log(w_j) − log(w_W), σ² )。两边的残差平方和除以 2σ²，就是似然贡献。

注意：这套似然和原始论文里用多项分布刻画离散判断的写法不完全一样，是工程上更常见的连续近似。好处是实现简单、跑得快；代价是它忽略了 1~9 整数刻度本身的离散性，权重极小时会有轻微偏差。做学术对比要按原论文的模型走，做业务赋权用这个够。

4.2 个体权重更新：Dirichlet 提议 + MH 接受率

主循环是一个 Metropolis-within-Gibbs：每一轮依次更新 K 个个体权重、一个全局权重、一个一致性参数，各步都用 MH，因为这三个条件分布都不共轭。

function [wNew, accepted] = update_wk(wOld, ab, aw, wStar, gamma, opts) %UPDATE_WK 单个决策人权重向量的 Metropolis 更新 % 先验：w^k ~ Dirichlet(gamma * wStar) % 似然：log(a_Bj) ~ N(log(w_B/w_j), sigma2) iB = find(ab == 1); % 最优准则下标 iW = find(aw == 1); % 最差准则下标 jB = setdiff(1:numel(ab), iB); % 与最优做比较的那些准则 jW = setdiff(1:numel(aw), iW); % 与最差做比较的那些准则 % Dirichlet 提议：均值落在当前点，kappa 决定步幅 alpha = opts.kappa * wOld; wNew = gamrnd(alpha, 1); wNew = wNew / sum(wNew); % 归一化回单纯形 % MH 接受率 = 后验比 × 提议密度比 logR = logpost(wNew, ab, aw, wStar, gamma, opts, iB, iW, jB, jW) ... - logpost(wOld, ab, aw, wStar, gamma, opts, iB, iW, jB, jW) ... + logq(wOld, wNew, opts.kappa) - logq(wNew, wOld, opts.kappa); accepted = log(rand) < logR; if ~accepted, wNew = wOld; end % 拒绝就原地不动 end function lp = logpost(w, ab, aw, wStar, gamma, opts, iB, iW, jB, jW) % 对数后验 = 对数先验 + 对数似然 lp = dirichlet_logpdf(w, gamma * wStar); % Dirichlet 先验 r1 = log(ab(jB)) - (log(w(iB)) - log(w(jB))); % 最优侧残差 r2 = log(aw(jW)) - (log(w(jW)) - log(w(iW))); % 最差侧残差 lp = lp - (sum(r1.^2) + sum(r2.^2)) / (2 * opts.sigma2); end function lq = logq(wTo, wFrom, kappa) % Dirichlet 提议密度 log q(wTo | wFrom) a = kappa * wFrom; lq = gammaln(kappa) - sum(gammaln(a)) + sum((a - 1) .* log(wTo)); end

几个点值得单独说。

alpha = opts.kappa * wOld是想让提议分布的中心刚好落在当前点上。κ 越大，Dirichlet 越尖锐，步幅越小；κ 越小，样本越散，步幅越大。经验值是让接受率落在 0.2~0.4 之间，对应的 κ 大概在 20~50。

gamrnd(alpha, 1)生成的是未归一化的 Gamma 样本，逐元素独立，再除以总和就落到单纯形上。这样提议出来的点自动满足 Σw = 1、w > 0，不需要额外加约束，也不用拒绝非法解。

logq那两项为什么必须留着？因为 Dirichlet 提议不对称——从 wOld 跳到 wNew 的概率，和从 wNew 跳回 wOld 的概率不是一回事。这是最容易写错的地方。我见过不少手写版本直接省掉这两项，链照样能跑、图照样能出，但收敛到的是一个错误的不变分布，权重会有系统性偏移，光看结果根本发现不了。

if ~accepted, wNew = wOld; end这行也别漏。拒绝的时候状态必须保留，写成wNew = wOld + 0或者干脆不赋值都会破坏细致平衡。

4.3 全局权重与一致性参数 γ 的更新

全局权重 w* 的先验取 Dirichlet(1)，也就是单纯形上的均匀分布；它的似然来自所有个体权重的先验密度连乘。展开之后

∏_k Dir(w^k | γ·w*) ∝ ∏_i w*_i^(γ·Σ_k w_i^k − K)

这个式子对 w* 不共轭，因为 Dirichlet 的归一化常数里也含 w*，所以还是得用 MH：

function [wNew, accepted] = update_wstar(wOld, W, gamma, opts) %UPDATE_WSTAR 全局权重 w* 的 Metropolis 更新 wNew = gamrnd(opts.kappa * wOld, 1); wNew = wNew / sum(wNew); lpo = dirichlet_logpdf(wOld, ones(size(wOld))) + sum_dir_like(wOld, W, gamma); lpn = dirichlet_logpdf(wNew, ones(size(wNew))) + sum_dir_like(wNew, W, gamma); logR = lpn - lpo + logq(wOld, wNew, opts.kappa) - logq(wNew, wOld, opts.kappa); accepted = log(rand) < logR; if ~accepted, wNew = wOld; end end function s = sum_dir_like(wStar, W, gamma) % 所有决策人先验密度对 w* 的贡献 s = 0; for k = 1:size(W, 2) s = s + dirichlet_logpdf(W(:,k), gamma * wStar); end end

γ 的更新用对数随机游走，因为 γ 必须为正：

function [gNew, accepted] = update_gamma(gOld, W, wStar, opts) %GAMMA 在 (1, gammaMax) 上做对数随机游走 gNew = exp(log(gOld) + 0.15 * randn); % 0.15 是步长 if gNew < 1 || gNew > opts.gammaMax % 越界直接拒绝 gNew = gOld; accepted = false; return; end lp = @(g) sum(arrayfun(@(k) dirichlet_logpdf(W(:,k), g * wStar), ... 1:size(W, 2))); logR = lp(gNew) - lp(gOld) + log(gNew) - log(gOld); % 末项为雅可比 accepted = log(rand) < logR; if ~accepted, gNew = gOld; end end

末项log(gNew) - log(gOld)是变换的雅可比。对数随机游走是在 log γ 空间里对称的，换回 γ 空间之后密度要乘上 1/γ，这一项不能省。同样是个"省了也能跑、跑出来是错的"的坑。

γ 的上界gammaMax是工程上的妥协。理论上 γ 可以取任意大，实际数据里 γ 一旦冲到几十，说明所有专家的判断几乎完全一致，这时候上界只是防止数值溢出。真出现 γ 贴着上界不动的情况，先别急着调大上界，回去看看数据是不是所有人填了同一份问卷。

4.4 收敛诊断：链没收敛时该看哪几个数

MCMC 跑完不能直接信输出，得先看三样东西。

接受率。bayes_bwm返回的 accRate 在 0.2~0.4 之间算健康。低于 0.1 说明步幅太大，链在原地打转，把 κ 调大；高于 0.6 说明步幅太小，一万轮走不了多远，把 κ 调小。

轨迹图。plot_results里画的是 γ 的迭代轨迹。健康的轨迹应该像一条快速抖动的毛线，没有明显的上升或下降趋势，也没有长时间贴在某一个值上。如果前 3000 轮一直在往下掉、之后才平稳，说明 burnIn 设短了，加到 5000 再试。

多链一致性。单链看不出问题，至少跑三条不同种子的链，比较各自 w* 的后验均值。三条链的差异应该在 0.01 量级以内。差异大就说明后验是多峰的，或者迭代次数远远不够。

想更严格一点，可以算 Gelman-Rubin 统计量 R̂：把每条链的样本分成前后两段，比较链间方差和链内方差。R̂ 小于 1.1 一般认为收敛。这份代码里没内置，但用summarize_posterior拿到的样本矩阵，自己写十行就能算出来，值不值得加看你的使用场景——做业务报告，轨迹图够用；发论文，R̂ 得有。

5. 避坑与排查：五类一跑就错的结果

5.1 权重全是 1/n，一点区分度都没有

现象：跑完一万轮，四个准则权重都接近 0.25，可信区间宽得能装下整个单纯形。

原因：AB 或 AW 里最优/最差那一位没填 1，find(ab == 1)返回空数组，jB就变成了全部准则，似然项里log(w(iB))直接报错或者算出 NaN，整个似然退化成常数，后验就只剩先验了。

解决：入口三行assert加上，特别是all(sum(AB==1,1)==1)这一条，能挡住 90% 的这类问题。另外检查是不是有专家把两个准则都填了 1。

5.2 接受率不到 5%，链几乎不动

现象：accRate打印出来 0.02，γ 的轨迹图是一条几乎水平的直线。

原因：sigma2设太小。似然项是-残差平方和 / (2*sigma2)，sigma2 = 0.01 的时候，稍微偏离一点的提议就会被压到 −Inf，全部被拒。

解决：先把 sigma2 放到 0.15 试跑，看接受率；再按 0.1、0.2、0.3 做一轮敏感性，看权重均值稳不稳。稳的话说明结论对噪声假设不敏感，报告里可以只写

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