简介面向多准则决策分析、贝叶斯统计建模与论文复现需求的开源 MATLAB 实现围绕贝叶斯最佳-最差方法Bayesian BWM展开通过概率分布刻画准则权重的不确定性可用于准则权重估计、方案排序评价以及相关实验的快速验证。压缩包共 8 个文件以 6 个 .m 脚本为主体辅以 1 个 markdown 说明和 1 个 txt 数据文本整体约 15KB体量轻便便于梳理模型调用、链收敛处理与结果绘图等环节。运行前需安装 JAGS 并正确配置环境变量示例中替换准则名称、最佳到其余准则比较向量等数据后即可开展自己的分析适合具备一定 MATLAB 基础的研究生、教师或算法实践者。已有 2174 人学习可作为课程作业、论文实验和决策分析原型搭建的参考帮助读者理解贝叶斯 BWM 从数据组织到结果解读的完整流程。1. 五个专家两套排序传统 BWM 藏起来的分歧有个做供应商评价的小团队五位专家各填一份最优最差问卷传统 BWM 跑出五套权重算术平均后排序是 ABC可其中两位专家的原始排序明明是 CAB——这种聚合等于把分歧直接抹平了。我第一次翻到 MATLAB 版贝叶斯 BWMBayesian 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 1w_j ≥ 0解完得到唯一的权重向量和一个一致性比率 ξ*。这套东西的好处是问的问题少2n−3 次比较比 AHP 的 n(n−1)/2 少一个量级坏处也很明显它只给你一个点估计。ξ* 是 0.08 还是 0.13你对这组权重的信任度是一样的没有区间、没有分布也没法回答准则 2 排第一的概率有多大这种问题。贝叶斯版本换了个思路把权重当成随机变量写成三层结构层次变量含义先验设定数据层AB、AW每位决策人的两组整数判断观测值不设先验个体层w^k第 k 位决策人的权重向量Dirichlet(γ·w*)群体层w*聚合后的全局权重Dirichlet(1)超参数γ决策人之间的一致性强度(1, γmax) 上均匀关键在于 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 Toolboxgamrnd可以用-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进去一路往下数据层main.m里把 AB、AW 两个 n×K 矩阵填好n 是准则数K 是决策人数模型层bayes_bwm.m接收矩阵初始化 w*、w^k、γ 三个状态量采样层循环里依次调update_wk→update_wstar→update_gamma每轮把状态推进一步展示层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三位专家都把质量选为最优、成本选为最差专家AB(C1)AB(C2)AB(C3)AB(C4)AW(C1)AW(C2)AW(C3)AW(C4)174131285265121394383141274读法专家 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(AB1,1)1), AB 每列必须有且仅有一个 1); assert(all(sum(AW1,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 后验均值 %.2f95%% 区间 [%.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]); endprctile(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*所以还是得用 MHfunction [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 里最优/最差那一位没填 1find(ab 1)返回空数组jB就变成了全部准则似然项里log(w(iB))直接报错或者算出 NaN整个似然退化成常数后验就只剩先验了。解决入口三行assert加上特别是all(sum(AB1,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 做一轮敏感性看权重均值稳不稳。稳的话说明结论对噪声假设不敏感报告里可以只写本文还有配套的精品资源点击获取