简介面向多准则决策、贝叶斯统计建模与MATLAB实现的研究生、教师及算法工程师这套开源代码围绕贝叶斯BWM提供可运行实现可用于指标权重推断、方案评价、论文复现及课程设计对多指标排序类问题具有参考价值。压缩包共8个文件约15KB以6个MATLAB源文件为主体辅以1个Markdown说明和1个TXT示例配置内容覆盖贝叶斯模型调用、结果绘图、MCMC链重构与收敛诊断等环节结构紧凑便于替换示例数据后开展自己的案例。环境上依赖JAGS完成概率采样读者可据此理解先验设定、后验计算与权重比较流程观察参数变化对结果的影响并积累贝叶斯实现的排错思路。目前已有2172人学习适合希望把BWM从理论公式落到可执行脚本、补齐实验代码与结果解读的读者参考。1. 贝叶斯BWM为什么“权重算出来了”不等于“排序站得住”把七个准则交给十位专家打分原始 BWM 十几分钟就能吐出一组权重和一张漂亮的排序表换一批专家排序就变了。多数人的第一反应是“专家不靠谱”真正的问题在于原始 BWM 只给一个点估计它告诉你 w10.31、w20.29却不告诉你这两个数的差距到底算不算差距。贝叶斯 BWM 把权重从待求参数改成随机变量用 Gamma 似然刻画专家打分的生成过程、用 Dirichlet 先验约束权重单纯形再靠 MCMC 把后验分布采出来排序于是从一张全序表变成一组带概率的 credal ranking。这套开源的 MATLAB 实现把模型、采样、置信度三段链路完整摊开改一个先验就能整条重跑。2. 原始 BWM 的压缩机制2n-3 次比较如何还原一组权重2.1 Best-to-Others 与 Others-to-Worst把 n(n-1)/2 个格子砍到两行AHP 式成对比较矩阵要求决策者填满上三角n 个准则就是 n(n-1)/2 次判断n8 时是 28 次专家填到后面基本靠惯性一致性比率随之崩掉。BWM 的做法是先在准则集里挑一个最好Best和一个最差Worst然后只问两类问题最好准则相对于其他每个准则重要多少倍其他每个准则相对于最差准则重要多少倍。这两个向量分别记作 a_B [a_B1 … a_Bn] 和 a_W [a_1W … a_nW]标度沿用 1 到 9 的整数。第一类有 n-1 个分量第二类也有 n-1 个其中 (Best, Worst) 这一对被问了两遍于是净比较次数是 2n-3。n8 时是 13 次比 28 次少一半还多这是 BWM 能被专家认真填完的根本原因。方法需要填的格子一致性处理输出形态全矩阵成对比较n(n-1)/2一致性比率 CR权重点估计原始 BWM2n-3一致性半径 ξ*权重点估计贝叶斯 BWM2n-3后验收缩 credal ranking权重联合分布两个向量的方向必须分清a_Bj 是“最好比 j 好多少”恒有 a_Bj ≥ 1a_jW 是“j 比最差好多少”同样恒有 a_jW ≥ 1。方向搞反是后面避坑章节里要专门讲的一类翻车。2.2 min-max 模型ξ 不是误差是一致性被允许的半径BWM 的理想关系是 a_Bj w_B / w_j 且 a_jW w_j / w_W。现实中专家做不到完全自洽于是把偏差收进一个统一的半径 ξmin ξ s.t. | w_B / w_j − a_Bj | ≤ ξ 对每个 j ≠ B | w_j / w_W − a_jW | ≤ ξ 对每个 j ≠ W Σ w_j 1, w_j 0注意 ξ 是作用在比值上的不是作用在权重上的所以它衡量的是一组权重能在多大程度上同时满足所有专家判断。ξ*0 意味着这套打分完全自洽ξ* 越大说明专家的判断越互相打架。工程上通常再算一个一致性比率 CR ξ*/CICI 是按标度查表得到的随机一致性指标CR 超过 0.1~0.2 就该回头找专家复谈。这里有个容易被忽略的点ξ 是拿“最坏的那一对”定义的它由单个最不一致的判断决定所以对个别异常打分极其敏感一个手滑填成 9 的分量就能把 ξ 顶上去。2.3 线性规划只能给出一个点单峰目标把不确定性抹平了上面的绝对值约束可以拆成两组线性不等式变量取 [w; ξ]直接交给 MATLAB 的 linprog。核心代码如下只保留了约束拼装部分function [w, xi] bwm_minmax(aB, aW) % aB: 1×n Best-to-OthersaW: 1×n Others-to-Worst n numel(aB); [~, iB] min(aB); % aBB 1 是最小值最小值位置即最佳准则 [~, iW] max(aW); % aWW 1 是最大值最大值位置即最差准则 A zeros(4*(n-1), n1); b zeros(4*(n-1), 1); r 0; for j 1:n if j ~ iB r r 1; A(r, [iB j n1]) [ 1, -aB(j), -1]; % wB - aBj*wj - xi 0 r r 1; A(r, [iB j n1]) [-1, aB(j), -1]; % -wB aBj*wj - xi 0 end if j ~ iW r r 1; A(r, [j iW n1]) [ 1, -aW(j), -1]; % wj - ajW*wW - xi 0 r r 1; A(r, [j iW n1]) [-1, aW(j), -1]; end end A A(1:r, :); b b(1:r); Aeq [ones(1,n), 0]; beq 1; % 权重和为 1xi 不参与该约束 lb [zeros(n,1); 0]; % 权重非负xi 非负 x linprog([zeros(n,1); 1], A, b, Aeq, beq, lb, [], ... optimoptions(linprog,Display,off)); w x(1:n); xi x(end); end参数说明输入两个长度均为 n 的行向量位置 iB 与 iW 处的元素必须是 1输出 w 是归一化权重xi 是一致性半径。每一对比较贡献两行不等式分别对应绝对值的正负两侧所以约束总数是 4(n-1)。目标函数系数只有 ξ 那一位是 1其余为 0因为权重本身不进入目标它们只是被约束住的自由变量。问题就出在这最后一句。线性规划的最优解是可行域的一个顶点它是一个确定的数没有任何“这组权重有多可信”的信息。哪怕把专家打分整体上下扰动一格linprog 也只是给你另一个顶点两个顶点之间没有可比的尺度。更麻烦的是 1 到 9 的标度本质上是序数性的语义刻度a_Bj7 和 a_Bj9 的差距并不等于 7 和 8 的差距把点估计的差值当成效果量去解读属于典型的过度解读。贝叶斯化要解决的就是这一层不是求一个最优点而是求权重在给定打分下的整个分布。3. 贝叶斯化把权重从“待求参数”改成“随机变量”之后3.1 似然函数a_Bj 为什么适合用 Gamma 分布承载原始 BWM 要求 a_Bj 精确等于 w_B/w_j贝叶斯版本改成“a_Bj 由 w_B/w_j 生成带随机扰动”。能承载这个角色的分布要满足三条取值必须为正比值不可能为负、右偏专家更容易高估而不是低估、并且期望或众数与 w_B/w_j 挂钩。Gamma 分布三条全中而且形状参数与尺度参数分开天然可以表达“这个专家有多稳”。常见的实现思路是把形状参数写成 κ·r、率参数写成 κ其中 r w_B/w_j 是模型期望的比值κ 控制离散程度a_Bj | w ~ Gamma(shape κ · w_B/w_j, rate κ) a_jW | w ~ Gamma(shape κ · w_j/w_W, rate κ)这样均值恰好是 w_B/w_j方差是 w_B/(κ·w_j)κ 越大分布越尖等于告诉模型“我更相信专家的数字”。这套参数化在文献里有若干变体有的把 κ 固定成常数有的把 κ 也当成超参数去采样思路是一致的形状参数里必须含有待估的权重比模型才有“从打分反推权重”的通道。权重的先验用 Dirichlet(α)α 取全 1 时就是单纯形上的均匀分布等价于“我对哪个准则更重要毫无先验”。这一步看着平淡实际很关键Dirichlet 先验自动把权重约束在单纯形上采样时不需要额外加 Σw1 的投影省掉了投影带来的边界偏差。3.2 先验与超参数 κ无信息不等于无影响很多人以为 α 全 1 就万事大吉真正左右后验形状的是 κ。κ 很大时似然极度尖锐后验几乎坍缩到 linprog 那个顶点附近贝叶斯化白做κ 很小时似然被压平后验几乎等于先验排序变成随机数。中间那段才是有效区间判断办法是做敏感性分析κ 取 0.1、1、10 各跑一遍看 credal ranking 的结论是否稳定。稳定的那部分结论才是数据支持的随 κ 翻脸的结论说明专家打分的信息量本来就不够。参数作用常用取值调错的症状αDirichlet 先验浓度全 1取大值会把权重往均匀拉κ似然离散度0.1 / 1 / 10 做敏感性太大后验坍缩太小后验发散nIter总迭代数1e4 ~ 1e5太短后验统计量不稳burnIn预热丢弃数总量的 20%~30%太短残留初值影响step随机游走步长burn-in 自适应接受率极端、自相关高3.3 后验没有闭式解共轭性断在哪一步把似然和先验乘起来取对数能得到一个可直接写进代码的目标函数。下面这段是权重空间下的对数后验逐项对应 Gamma 密度的对数function lp log_post(w, aB, aW, iB, iW, kappa, alpha) % w: 归一化权重列向量kappa: 似然离散度alpha: Dirichlet 浓度向量 if any(w 0), lp -inf; return; end lp 0; for j 1:numel(w) if j ~ iB % Best-to-Others 这一支 r w(iB) / w(j); k kappa * r; % Gamma 形状参数 lp lp k*log(kappa) (k-1)*log(aB(j)) - kappa*aB(j) - gammaln(k); end if j ~ iW % Others-to-Worst 这一支 r w(j) / w(iW); k kappa * r; lp lp k*log(kappa) (k-1)*log(aW(j)) - kappa*aW(j) - gammaln(k); end end lp lp sum((alpha - 1) .* log(w)); % Dirichlet 先验 end三个细节值得盯住。第一Gamma 密度的归一化项是 gammaln(k)而 k 里含 w所以这一项不能当常数丢掉丢掉会让后验系统性偏斜这是手写贝叶斯代码最常见的隐形 bug。第二j iB 和 j iW 必须跳过因为 a_BB1、a_WW1 是定义而不是观测塞进似然等于给权重加了一个不该有的锚。第三循环里两支的结构不对称Best 那一支的分子是固定的 w_BWorst 那一支的分母是固定的 w_W写反了整条链会收敛到一个错误分布而且看上去还挺像回事。共轭性就断在 k 依赖 w 这一点上。把某个 w_j 单独拎出来看条件后验里会同时出现 exp(A/w_j)、exp(B·w_j)、gammaln(c/w_j)、gammaln(c·w_j) 四类项凑不成任何标准分布所以 Gibbs 采样没法直接“抽一发”必须混进 Metropolis-Hastings 或者改用 HMC。这就是为什么这套实现不能像 Dirichlet-multinomial 那样几行写完底层模型的非共轭性决定了采样器必须自己搭。4. 采样与 Credal Ranking置信度怎么从一条马尔可夫链里算出来4.1 softmax 变换与 Metropolis-within-Gibbs 的每一步权重被约束在单纯形上直接在单纯形上做随机游走会很别扭你对某个分量加一点别的分量就得按比例减接受率计算还要额外处理。更干净的做法是换到无约束空间用 softmax 把任意实数向量 z 映射成权重 w。代价是必须补一个雅可比项——softmax 的雅可比行列式是 ∏ w_i所以无约束空间下的对数目标密度是 log_posterior(w) Σ log w_i。很多人以为 Dirichlet(1) 先验已经把“均匀”表达完了于是漏掉这一项链照样收敛、图照样好看但后验是偏的这就是后面要讲的第一号坑。Metropolis-within-Gibbs 是折中方案外层按坐标轮转内层对单个 z_j 做对称的 Gaussian 随机游走。对称提议的接受率简化为目标密度之比不用再算提议分布的比值实现上只要两行function lz log_target(z, aB, aW, iB, iW, kappa, alpha) z z - max(z); % 数值稳定防止 exp 溢出 w exp(z) / sum(exp(z)); % softmax 映射到单纯形 lz log_post(w, aB, aW, iB, iW, kappa, alpha) sum(log(w)); % 雅可比项 end主循环里有两处设计取舍值得说明。一是步长只在 burn-in 阶段自适应进入采样阶段后必须冻结边采样边改步长会破坏细致平衡链收敛到的就不是目标分布了这种错误极难从 trace 图上看出来。二是多链不同初值靠链间差异判断收敛单链跑得再久也证明不了什么。nIter 30000; burnIn 9000; z zeros(1, n); step 0.1 * ones(1, n); lz log_target(z, aB, aW, iB, iW, kappa, alpha); Ws zeros(nIter - burnIn, n); acc zeros(1,n); tries zeros(1,n); for it 1:nIter for j 1:n zc z; zc(j) z(j) step(j) * randn; % 对称随机游走 lzc log_target(zc, aB, aW, iB, iW, kappa, alpha); if log(rand) lzc - lz z zc; lz lzc; acc(j) acc(j) 1; end tries(j) tries(j) 1; end if it burnIn mod(it, 200) 0 % 只在预热期自适应 rate acc ./ tries; step(rate 0.2) step(rate 0.2) * 0.8; step(rate 0.5) step(rate 0.5) * 1.2; acc(:) 0; tries(:) 0; end if it burnIn w exp(z - max(z)); Ws(it - burnIn, :) w / sum(w); end end每个坐标的目标接受率定在 0.3 附近比较合适低于 0.2 就缩步长、高于 0.5 就放步长这种粗调在维度不高时够用。n 个准则、三万次迭代MATLAB 里通常跑几秒到几十秒比接一个外部概率编程环境省事得多。4.2 后验均值、可信区间与 ξ 的后验分布采完之后权重不再是一个向量而是一个样本矩阵 Ws样本数 × n。逐列取均值得到后验均值权重逐列取 2.5% 和 97.5% 分位数得到 95% 可信区间。区间宽窄本身就是信息某个准则的区间跨了 0.1说明专家在这个准则上的意见本来就没收敛再怎么算也排不出确定名次。一致性半径也能拿到分布对每一行样本重新算一次 max|w_B/w_j − a_Bj| 与 max|w_j/w_W − a_jW|取两者的较大值就得到 ξ 的一条后验样本。这样得到的是 ξ 的完整分布而不是一个数可以直接回答“这组打分不大可能一致的概率是多少”。相比之下原始 BWM 只给一个 ξ*你没有任何尺度判断它偏大还是偏小只能死记 CR 0.1 这个阈值。4.3 credal ranking把全序拆成“有把握”和“不可比”有样本矩阵之后任意两个准则谁更重要的概率就是一个计数问题P(w_i w_j) 等于样本中第 i 列大于第 j 列的频率。这个概率矩阵就是 credal ranking 的全部原料剩下的只是阈值选择。function [P, pairs] credal_ranking(Ws, thr) % Ws: nSample×n 的后验样本矩阵thr: 认定“优于”的概率阈值 n size(Ws, 2); P zeros(n); pairs {}; for i 1:n for j 1:n if i ~ j P(i,j) mean(Ws(:,i) Ws(:,j)); % 后验概率 end end end for i 1:n for j i1:n if P(i,j) thr pairs{end1} sprintf(w%d w%d (%.2f), i, j, P(i,j)); elseif P(j,i) thr pairs{end1} sprintf(w%d w%d (%.2f), j, i, P(j,i)); end end end end阈值取 0.5 时几乎每对准则都能分出胜负结果接近一个全序但代价是把 51% 对 49% 这种毫无意义的比较也当成了结论。把阈值抬到 0.75 或 0.9一部分准则对就变成“不可比”输出的是一张偏序表。这不是缺陷而是特性两个准则的后验分布大面积重叠时硬排出一个名次只是在制造虚假精度。如果阈值 0.5 下出现了 AB、BC、CA 的环说明这三个准则在数据上根本分不开该做的是回去检查打分或者合并准则而不是换个随机种子重跑。5. 避坑与排查跑贝叶斯 BWM 最容易翻车的五个地方5.1 漏掉 softmax 雅可比链收敛了但后验是偏的现象trace 图平稳接受率正常但后验均值权重与 linprog 的点估计系统性偏离且偏离方向随 κ 变化。原因换到无约束空间采样时只用了 log_post(w)忘了加 ∏ w_i 这一项雅可比。Dirichlet(1) 先验在权重空间是均匀的映射到 z 空间就变成 ∏ w_i 的密度漏掉它等于偷偷换了一个更偏向均匀权重的先验。解决把雅可比项写进 log_target并且做一次定量校验——固定 κ 到一个很大的值比如 100后验均值应当收敛到 linprog 的解附近对不上就说明目标密度写错了。5.2 把 a_BB 和 a_WW 当观测塞进似然现象一致性半径 ξ 的后验几乎全为 0权重区间窄得不像话结论好得可疑。原因a_BB1 和 a_WW1 是定义值不是专家给出的判断。把它们当数据放进似然相当于额外加了两条“w_B/w_B 必须等于 1”的锚约束人为压低了后验方差。解决循环里显式跳过 j iB 和 j iW更稳的做法是在进入采样器之前就把这两个位置从向量里剔掉让数据结构上不可能出错。5.3 步长写死接受率掉到个位数现象样本自相关爆表有效样本量可能只有几百credal 概率在 0.5 附近来回抖动换个随机种子结论就变。原因n 个准则的尺度不同w 小的分量在 z 空间需要的步长也小统一写死一个 step 会让某些坐标几乎永远被拒绝。解决每个坐标一个步长burn-in 阶段按接受率自适应采样阶段冻结同时开 3~4 条不同初值的链比较链间方差。5.4 aW 方向填反最差准则反而权重最高现象后验均值里本该最差的准则排到了第一而且 credal ranking 置信度还很“坚定”。原因a_Bj 是“最好比 j 好多少”a_jW 是“j 比最差好多少”两者都应当 ≥1但方向相反。如果按 a_W 的形式去填 a_jW 的倒数模型看到的是“j 比最差差多少”整个比例关系被翻转。解决在数据入口处做一次断言检查所有 a_B 分量 ≥1、所有 a_W 分量 ≥1、且 a_B(iW) a_W(iB)。这三条过不了就不许进采样器。5.5 把 credal 概率当置信水平报出去现象报告里写“准则 A 优于准则 B 的概率为 95%”被追问这是什么意义上的概率时答不上来。原因P(w_A w_B)0.95 是给定数据和模型后的后验概率不是重复抽样意义下的频率置信度而且它依赖 κ 的取值和先验设定换个 κ 就可能变成 0.7。解决报整条后验分布和区间附一张 κ 的敏感性表结论只写“在 κ1 的模型设定下A 优于 B 的后验概率为 0.95”把模型假设显式写在结论旁边。6. 设计取舍手写采样器、Gamma 似然与“慢一点换透明”这套实现最容易被质疑的地方是为什么不直接接到现成的概率编程环境里用 HMC 或 NUTS 采样几十行就能搞定还自带诊断答案在于可读性和可改性之间的权衡。HMC 的采样效率在手写 MH 的十倍以上但它把模型编译成一个黑匣子用户想改先验、想换似然、想看某个中间量都得重新学一套建模语言而 MATLAB 里这套 Metropolis-within-Gibbs 的每一步都摊在明面上改 κ、换 Dirichlet 参数、加一个新的似然分支改动量就是几行循环。对教学和小规模准则集n 通常不超过 15来说几十秒的运行时间换来的透明是划算的。这是典型的“用算力换可解释”不是性能最优解。第二个取舍在 Gamma 似然的选择上。用正态分布去拟合 log 比值也能跑而且共轭性更好写但正态在 log 尺度上对称意味着模型认为专家高估和低估的概率一样大这与实际打分行为不符——面对 1 到 9 的整数标度专家更容易往大了报。Gamma 的正偏形态更贴近真实数据代价是形状参数含 w 导致的非共轭性必须自己写采样器。换句话说5.1 那条雅可比坑、5.3 那条步长坑本质都是为“似然更贴近现实”付出的工程税。第三个取舍是 credal ranking本文还有配套的精品资源点击获取 
长河君 · 2026/10/29 15:40:31 · 阅读 1.2k 贝叶斯BWM：权重不确定性、credal ranking与MATLAB实现避坑
简介面向多准则决策、贝叶斯统计建模与MATLAB实现的研究生、教师及算法工程师这套开源代码围绕贝叶斯BWM提供可运行实现可用于指标权重推断、方案评价、论文复现及课程设计对多指标排序类问题具有参考价值。压缩包共8个文件约15KB以6个MATLAB源文件为主体辅以1个Markdown说明和1个TXT示例配置内容覆盖贝叶斯模型调用、结果绘图、MCMC链重构与收敛诊断等环节结构紧凑便于替换示例数据后开展自己的案例。环境上依赖JAGS完成概率采样读者可据此理解先验设定、后验计算与权重比较流程观察参数变化对结果的影响并积累贝叶斯实现的排错思路。目前已有2172人学习适合希望把BWM从理论公式落到可执行脚本、补齐实验代码与结果解读的读者参考。1. 贝叶斯BWM为什么“权重算出来了”不等于“排序站得住”把七个准则交给十位专家打分原始 BWM 十几分钟就能吐出一组权重和一张漂亮的排序表换一批专家排序就变了。多数人的第一反应是“专家不靠谱”真正的问题在于原始 BWM 只给一个点估计它告诉你 w10.31、w20.29却不告诉你这两个数的差距到底算不算差距。贝叶斯 BWM 把权重从待求参数改成随机变量用 Gamma 似然刻画专家打分的生成过程、用 Dirichlet 先验约束权重单纯形再靠 MCMC 把后验分布采出来排序于是从一张全序表变成一组带概率的 credal ranking。这套开源的 MATLAB 实现把模型、采样、置信度三段链路完整摊开改一个先验就能整条重跑。2. 原始 BWM 的压缩机制2n-3 次比较如何还原一组权重2.1 Best-to-Others 与 Others-to-Worst把 n(n-1)/2 个格子砍到两行AHP 式成对比较矩阵要求决策者填满上三角n 个准则就是 n(n-1)/2 次判断n8 时是 28 次专家填到后面基本靠惯性一致性比率随之崩掉。BWM 的做法是先在准则集里挑一个最好Best和一个最差Worst然后只问两类问题最好准则相对于其他每个准则重要多少倍其他每个准则相对于最差准则重要多少倍。这两个向量分别记作 a_B [a_B1 … a_Bn] 和 a_W [a_1W … a_nW]标度沿用 1 到 9 的整数。第一类有 n-1 个分量第二类也有 n-1 个其中 (Best, Worst) 这一对被问了两遍于是净比较次数是 2n-3。n8 时是 13 次比 28 次少一半还多这是 BWM 能被专家认真填完的根本原因。方法需要填的格子一致性处理输出形态全矩阵成对比较n(n-1)/2一致性比率 CR权重点估计原始 BWM2n-3一致性半径 ξ*权重点估计贝叶斯 BWM2n-3后验收缩 credal ranking权重联合分布两个向量的方向必须分清a_Bj 是“最好比 j 好多少”恒有 a_Bj ≥ 1a_jW 是“j 比最差好多少”同样恒有 a_jW ≥ 1。方向搞反是后面避坑章节里要专门讲的一类翻车。2.2 min-max 模型ξ 不是误差是一致性被允许的半径BWM 的理想关系是 a_Bj w_B / w_j 且 a_jW w_j / w_W。现实中专家做不到完全自洽于是把偏差收进一个统一的半径 ξmin ξ s.t. | w_B / w_j − a_Bj | ≤ ξ 对每个 j ≠ B | w_j / w_W − a_jW | ≤ ξ 对每个 j ≠ W Σ w_j 1, w_j 0注意 ξ 是作用在比值上的不是作用在权重上的所以它衡量的是一组权重能在多大程度上同时满足所有专家判断。ξ*0 意味着这套打分完全自洽ξ* 越大说明专家的判断越互相打架。工程上通常再算一个一致性比率 CR ξ*/CICI 是按标度查表得到的随机一致性指标CR 超过 0.1~0.2 就该回头找专家复谈。这里有个容易被忽略的点ξ 是拿“最坏的那一对”定义的它由单个最不一致的判断决定所以对个别异常打分极其敏感一个手滑填成 9 的分量就能把 ξ 顶上去。2.3 线性规划只能给出一个点单峰目标把不确定性抹平了上面的绝对值约束可以拆成两组线性不等式变量取 [w; ξ]直接交给 MATLAB 的 linprog。核心代码如下只保留了约束拼装部分function [w, xi] bwm_minmax(aB, aW) % aB: 1×n Best-to-OthersaW: 1×n Others-to-Worst n numel(aB); [~, iB] min(aB); % aBB 1 是最小值最小值位置即最佳准则 [~, iW] max(aW); % aWW 1 是最大值最大值位置即最差准则 A zeros(4*(n-1), n1); b zeros(4*(n-1), 1); r 0; for j 1:n if j ~ iB r r 1; A(r, [iB j n1]) [ 1, -aB(j), -1]; % wB - aBj*wj - xi 0 r r 1; A(r, [iB j n1]) [-1, aB(j), -1]; % -wB aBj*wj - xi 0 end if j ~ iW r r 1; A(r, [j iW n1]) [ 1, -aW(j), -1]; % wj - ajW*wW - xi 0 r r 1; A(r, [j iW n1]) [-1, aW(j), -1]; end end A A(1:r, :); b b(1:r); Aeq [ones(1,n), 0]; beq 1; % 权重和为 1xi 不参与该约束 lb [zeros(n,1); 0]; % 权重非负xi 非负 x linprog([zeros(n,1); 1], A, b, Aeq, beq, lb, [], ... optimoptions(linprog,Display,off)); w x(1:n); xi x(end); end参数说明输入两个长度均为 n 的行向量位置 iB 与 iW 处的元素必须是 1输出 w 是归一化权重xi 是一致性半径。每一对比较贡献两行不等式分别对应绝对值的正负两侧所以约束总数是 4(n-1)。目标函数系数只有 ξ 那一位是 1其余为 0因为权重本身不进入目标它们只是被约束住的自由变量。问题就出在这最后一句。线性规划的最优解是可行域的一个顶点它是一个确定的数没有任何“这组权重有多可信”的信息。哪怕把专家打分整体上下扰动一格linprog 也只是给你另一个顶点两个顶点之间没有可比的尺度。更麻烦的是 1 到 9 的标度本质上是序数性的语义刻度a_Bj7 和 a_Bj9 的差距并不等于 7 和 8 的差距把点估计的差值当成效果量去解读属于典型的过度解读。贝叶斯化要解决的就是这一层不是求一个最优点而是求权重在给定打分下的整个分布。3. 贝叶斯化把权重从“待求参数”改成“随机变量”之后3.1 似然函数a_Bj 为什么适合用 Gamma 分布承载原始 BWM 要求 a_Bj 精确等于 w_B/w_j贝叶斯版本改成“a_Bj 由 w_B/w_j 生成带随机扰动”。能承载这个角色的分布要满足三条取值必须为正比值不可能为负、右偏专家更容易高估而不是低估、并且期望或众数与 w_B/w_j 挂钩。Gamma 分布三条全中而且形状参数与尺度参数分开天然可以表达“这个专家有多稳”。常见的实现思路是把形状参数写成 κ·r、率参数写成 κ其中 r w_B/w_j 是模型期望的比值κ 控制离散程度a_Bj | w ~ Gamma(shape κ · w_B/w_j, rate κ) a_jW | w ~ Gamma(shape κ · w_j/w_W, rate κ)这样均值恰好是 w_B/w_j方差是 w_B/(κ·w_j)κ 越大分布越尖等于告诉模型“我更相信专家的数字”。这套参数化在文献里有若干变体有的把 κ 固定成常数有的把 κ 也当成超参数去采样思路是一致的形状参数里必须含有待估的权重比模型才有“从打分反推权重”的通道。权重的先验用 Dirichlet(α)α 取全 1 时就是单纯形上的均匀分布等价于“我对哪个准则更重要毫无先验”。这一步看着平淡实际很关键Dirichlet 先验自动把权重约束在单纯形上采样时不需要额外加 Σw1 的投影省掉了投影带来的边界偏差。3.2 先验与超参数 κ无信息不等于无影响很多人以为 α 全 1 就万事大吉真正左右后验形状的是 κ。κ 很大时似然极度尖锐后验几乎坍缩到 linprog 那个顶点附近贝叶斯化白做κ 很小时似然被压平后验几乎等于先验排序变成随机数。中间那段才是有效区间判断办法是做敏感性分析κ 取 0.1、1、10 各跑一遍看 credal ranking 的结论是否稳定。稳定的那部分结论才是数据支持的随 κ 翻脸的结论说明专家打分的信息量本来就不够。参数作用常用取值调错的症状αDirichlet 先验浓度全 1取大值会把权重往均匀拉κ似然离散度0.1 / 1 / 10 做敏感性太大后验坍缩太小后验发散nIter总迭代数1e4 ~ 1e5太短后验统计量不稳burnIn预热丢弃数总量的 20%~30%太短残留初值影响step随机游走步长burn-in 自适应接受率极端、自相关高3.3 后验没有闭式解共轭性断在哪一步把似然和先验乘起来取对数能得到一个可直接写进代码的目标函数。下面这段是权重空间下的对数后验逐项对应 Gamma 密度的对数function lp log_post(w, aB, aW, iB, iW, kappa, alpha) % w: 归一化权重列向量kappa: 似然离散度alpha: Dirichlet 浓度向量 if any(w 0), lp -inf; return; end lp 0; for j 1:numel(w) if j ~ iB % Best-to-Others 这一支 r w(iB) / w(j); k kappa * r; % Gamma 形状参数 lp lp k*log(kappa) (k-1)*log(aB(j)) - kappa*aB(j) - gammaln(k); end if j ~ iW % Others-to-Worst 这一支 r w(j) / w(iW); k kappa * r; lp lp k*log(kappa) (k-1)*log(aW(j)) - kappa*aW(j) - gammaln(k); end end lp lp sum((alpha - 1) .* log(w)); % Dirichlet 先验 end三个细节值得盯住。第一Gamma 密度的归一化项是 gammaln(k)而 k 里含 w所以这一项不能当常数丢掉丢掉会让后验系统性偏斜这是手写贝叶斯代码最常见的隐形 bug。第二j iB 和 j iW 必须跳过因为 a_BB1、a_WW1 是定义而不是观测塞进似然等于给权重加了一个不该有的锚。第三循环里两支的结构不对称Best 那一支的分子是固定的 w_BWorst 那一支的分母是固定的 w_W写反了整条链会收敛到一个错误分布而且看上去还挺像回事。共轭性就断在 k 依赖 w 这一点上。把某个 w_j 单独拎出来看条件后验里会同时出现 exp(A/w_j)、exp(B·w_j)、gammaln(c/w_j)、gammaln(c·w_j) 四类项凑不成任何标准分布所以 Gibbs 采样没法直接“抽一发”必须混进 Metropolis-Hastings 或者改用 HMC。这就是为什么这套实现不能像 Dirichlet-multinomial 那样几行写完底层模型的非共轭性决定了采样器必须自己搭。4. 采样与 Credal Ranking置信度怎么从一条马尔可夫链里算出来4.1 softmax 变换与 Metropolis-within-Gibbs 的每一步权重被约束在单纯形上直接在单纯形上做随机游走会很别扭你对某个分量加一点别的分量就得按比例减接受率计算还要额外处理。更干净的做法是换到无约束空间用 softmax 把任意实数向量 z 映射成权重 w。代价是必须补一个雅可比项——softmax 的雅可比行列式是 ∏ w_i所以无约束空间下的对数目标密度是 log_posterior(w) Σ log w_i。很多人以为 Dirichlet(1) 先验已经把“均匀”表达完了于是漏掉这一项链照样收敛、图照样好看但后验是偏的这就是后面要讲的第一号坑。Metropolis-within-Gibbs 是折中方案外层按坐标轮转内层对单个 z_j 做对称的 Gaussian 随机游走。对称提议的接受率简化为目标密度之比不用再算提议分布的比值实现上只要两行function lz log_target(z, aB, aW, iB, iW, kappa, alpha) z z - max(z); % 数值稳定防止 exp 溢出 w exp(z) / sum(exp(z)); % softmax 映射到单纯形 lz log_post(w, aB, aW, iB, iW, kappa, alpha) sum(log(w)); % 雅可比项 end主循环里有两处设计取舍值得说明。一是步长只在 burn-in 阶段自适应进入采样阶段后必须冻结边采样边改步长会破坏细致平衡链收敛到的就不是目标分布了这种错误极难从 trace 图上看出来。二是多链不同初值靠链间差异判断收敛单链跑得再久也证明不了什么。nIter 30000; burnIn 9000; z zeros(1, n); step 0.1 * ones(1, n); lz log_target(z, aB, aW, iB, iW, kappa, alpha); Ws zeros(nIter - burnIn, n); acc zeros(1,n); tries zeros(1,n); for it 1:nIter for j 1:n zc z; zc(j) z(j) step(j) * randn; % 对称随机游走 lzc log_target(zc, aB, aW, iB, iW, kappa, alpha); if log(rand) lzc - lz z zc; lz lzc; acc(j) acc(j) 1; end tries(j) tries(j) 1; end if it burnIn mod(it, 200) 0 % 只在预热期自适应 rate acc ./ tries; step(rate 0.2) step(rate 0.2) * 0.8; step(rate 0.5) step(rate 0.5) * 1.2; acc(:) 0; tries(:) 0; end if it burnIn w exp(z - max(z)); Ws(it - burnIn, :) w / sum(w); end end每个坐标的目标接受率定在 0.3 附近比较合适低于 0.2 就缩步长、高于 0.5 就放步长这种粗调在维度不高时够用。n 个准则、三万次迭代MATLAB 里通常跑几秒到几十秒比接一个外部概率编程环境省事得多。4.2 后验均值、可信区间与 ξ 的后验分布采完之后权重不再是一个向量而是一个样本矩阵 Ws样本数 × n。逐列取均值得到后验均值权重逐列取 2.5% 和 97.5% 分位数得到 95% 可信区间。区间宽窄本身就是信息某个准则的区间跨了 0.1说明专家在这个准则上的意见本来就没收敛再怎么算也排不出确定名次。一致性半径也能拿到分布对每一行样本重新算一次 max|w_B/w_j − a_Bj| 与 max|w_j/w_W − a_jW|取两者的较大值就得到 ξ 的一条后验样本。这样得到的是 ξ 的完整分布而不是一个数可以直接回答“这组打分不大可能一致的概率是多少”。相比之下原始 BWM 只给一个 ξ*你没有任何尺度判断它偏大还是偏小只能死记 CR 0.1 这个阈值。4.3 credal ranking把全序拆成“有把握”和“不可比”有样本矩阵之后任意两个准则谁更重要的概率就是一个计数问题P(w_i w_j) 等于样本中第 i 列大于第 j 列的频率。这个概率矩阵就是 credal ranking 的全部原料剩下的只是阈值选择。function [P, pairs] credal_ranking(Ws, thr) % Ws: nSample×n 的后验样本矩阵thr: 认定“优于”的概率阈值 n size(Ws, 2); P zeros(n); pairs {}; for i 1:n for j 1:n if i ~ j P(i,j) mean(Ws(:,i) Ws(:,j)); % 后验概率 end end end for i 1:n for j i1:n if P(i,j) thr pairs{end1} sprintf(w%d w%d (%.2f), i, j, P(i,j)); elseif P(j,i) thr pairs{end1} sprintf(w%d w%d (%.2f), j, i, P(j,i)); end end end end阈值取 0.5 时几乎每对准则都能分出胜负结果接近一个全序但代价是把 51% 对 49% 这种毫无意义的比较也当成了结论。把阈值抬到 0.75 或 0.9一部分准则对就变成“不可比”输出的是一张偏序表。这不是缺陷而是特性两个准则的后验分布大面积重叠时硬排出一个名次只是在制造虚假精度。如果阈值 0.5 下出现了 AB、BC、CA 的环说明这三个准则在数据上根本分不开该做的是回去检查打分或者合并准则而不是换个随机种子重跑。5. 避坑与排查跑贝叶斯 BWM 最容易翻车的五个地方5.1 漏掉 softmax 雅可比链收敛了但后验是偏的现象trace 图平稳接受率正常但后验均值权重与 linprog 的点估计系统性偏离且偏离方向随 κ 变化。原因换到无约束空间采样时只用了 log_post(w)忘了加 ∏ w_i 这一项雅可比。Dirichlet(1) 先验在权重空间是均匀的映射到 z 空间就变成 ∏ w_i 的密度漏掉它等于偷偷换了一个更偏向均匀权重的先验。解决把雅可比项写进 log_target并且做一次定量校验——固定 κ 到一个很大的值比如 100后验均值应当收敛到 linprog 的解附近对不上就说明目标密度写错了。5.2 把 a_BB 和 a_WW 当观测塞进似然现象一致性半径 ξ 的后验几乎全为 0权重区间窄得不像话结论好得可疑。原因a_BB1 和 a_WW1 是定义值不是专家给出的判断。把它们当数据放进似然相当于额外加了两条“w_B/w_B 必须等于 1”的锚约束人为压低了后验方差。解决循环里显式跳过 j iB 和 j iW更稳的做法是在进入采样器之前就把这两个位置从向量里剔掉让数据结构上不可能出错。5.3 步长写死接受率掉到个位数现象样本自相关爆表有效样本量可能只有几百credal 概率在 0.5 附近来回抖动换个随机种子结论就变。原因n 个准则的尺度不同w 小的分量在 z 空间需要的步长也小统一写死一个 step 会让某些坐标几乎永远被拒绝。解决每个坐标一个步长burn-in 阶段按接受率自适应采样阶段冻结同时开 3~4 条不同初值的链比较链间方差。5.4 aW 方向填反最差准则反而权重最高现象后验均值里本该最差的准则排到了第一而且 credal ranking 置信度还很“坚定”。原因a_Bj 是“最好比 j 好多少”a_jW 是“j 比最差好多少”两者都应当 ≥1但方向相反。如果按 a_W 的形式去填 a_jW 的倒数模型看到的是“j 比最差差多少”整个比例关系被翻转。解决在数据入口处做一次断言检查所有 a_B 分量 ≥1、所有 a_W 分量 ≥1、且 a_B(iW) a_W(iB)。这三条过不了就不许进采样器。5.5 把 credal 概率当置信水平报出去现象报告里写“准则 A 优于准则 B 的概率为 95%”被追问这是什么意义上的概率时答不上来。原因P(w_A w_B)0.95 是给定数据和模型后的后验概率不是重复抽样意义下的频率置信度而且它依赖 κ 的取值和先验设定换个 κ 就可能变成 0.7。解决报整条后验分布和区间附一张 κ 的敏感性表结论只写“在 κ1 的模型设定下A 优于 B 的后验概率为 0.95”把模型假设显式写在结论旁边。6. 设计取舍手写采样器、Gamma 似然与“慢一点换透明”这套实现最容易被质疑的地方是为什么不直接接到现成的概率编程环境里用 HMC 或 NUTS 采样几十行就能搞定还自带诊断答案在于可读性和可改性之间的权衡。HMC 的采样效率在手写 MH 的十倍以上但它把模型编译成一个黑匣子用户想改先验、想换似然、想看某个中间量都得重新学一套建模语言而 MATLAB 里这套 Metropolis-within-Gibbs 的每一步都摊在明面上改 κ、换 Dirichlet 参数、加一个新的似然分支改动量就是几行循环。对教学和小规模准则集n 通常不超过 15来说几十秒的运行时间换来的透明是划算的。这是典型的“用算力换可解释”不是性能最优解。第二个取舍在 Gamma 似然的选择上。用正态分布去拟合 log 比值也能跑而且共轭性更好写但正态在 log 尺度上对称意味着模型认为专家高估和低估的概率一样大这与实际打分行为不符——面对 1 到 9 的整数标度专家更容易往大了报。Gamma 的正偏形态更贴近真实数据代价是形状参数含 w 导致的非共轭性必须自己写采样器。换句话说5.1 那条雅可比坑、5.3 那条步长坑本质都是为“似然更贴近现实”付出的工程税。第三个取舍是 credal ranking本文还有配套的精品资源点击获取