简介：面向多准则决策、贝叶斯统计建模与MATLAB实现的研究生、教师及算法工程师，这套开源代码围绕贝叶斯BWM提供可运行实现，可用于指标权重推断、方案评价、论文复现及课程设计，对多指标排序类问题具有参考价值。压缩包共8个文件，约15KB，以6个MATLAB源文件为主体，辅以1个Markdown说明和1个TXT示例配置；内容覆盖贝叶斯模型调用、结果绘图、MCMC链重构与收敛诊断等环节，结构紧凑，便于替换示例数据后开展自己的案例。环境上依赖JAGS完成概率采样，读者可据此理解先验设定、后验计算与权重比较流程，观察参数变化对结果的影响，并积累贝叶斯实现的排错思路。目前已有2172人学习，适合希望把BWM从理论公式落到可执行脚本、补齐实验代码与结果解读的读者参考。

1. 贝叶斯BWM：为什么“权重算出来了”不等于“排序站得住”

把七个准则交给十位专家打分，原始 BWM 十几分钟就能吐出一组权重和一张漂亮的排序表；换一批专家，排序就变了。多数人的第一反应是“专家不靠谱”，真正的问题在于原始 BWM 只给一个点估计：它告诉你 w1=0.31、w2=0.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 次判断，n=8 时是 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。n=8 时是 13 次，比 28 次少一半还多，这是 BWM 能被专家认真填完的根本原因。

两个向量的方向必须分清：a_Bj 是“最好比 j 好多少”，恒有 a_Bj ≥ 1；a_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 = ξ*/CI，CI 是按标度查表得到的随机一致性指标，CR 超过 0.1~0.2 就该回头找专家复谈。这里有个容易被忽略的点：ξ 是拿“最坏的那一对”定义的，它由单个最不一致的判断决定，所以对个别异常打分极其敏感，一个手滑填成 9 的分量就能把 ξ 顶上去。

2.3 线性规划只能给出一个点：单峰目标把不确定性抹平了

上面的绝对值约束可以拆成两组线性不等式，变量取 [w; ξ]，直接交给 MATLAB 的 linprog。核心代码如下，只保留了约束拼装部分：

function [w, xi] = bwm_minmax(aB, aW) % aB: 1×n Best-to-Others；aW: 1×n Others-to-Worst n = numel(aB); [~, iB] = min(aB); % aBB = 1 是最小值，最小值位置即最佳准则 [~, iW] = max(aW); % aWW = 1 是最大值，最大值位置即最差准则 A = zeros(4*(n-1), n+1); b = zeros(4*(n-1), 1); r = 0; for j = 1:n if j ~= iB r = r + 1; A(r, [iB j n+1]) = [ 1, -aB(j), -1]; % wB - aBj*wj - xi <= 0 r = r + 1; A(r, [iB j n+1]) = [-1, aB(j), -1]; % -wB + aBj*wj - xi <= 0 end if j ~= iW r = r + 1; A(r, [j iW n+1]) = [ 1, -aW(j), -1]; % wj - ajW*wW - xi <= 0 r = r + 1; A(r, [j iW n+1]) = [-1, aW(j), -1]; end end A = A(1:r, :); b = b(1:r); Aeq = [ones(1,n), 0]; beq = 1; % 权重和为 1，xi 不参与该约束 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_Bj=7 和 a_Bj=9 的差距并不等于 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 先验自动把权重约束在单纯形上，采样时不需要额外加 Σw=1 的投影，省掉了投影带来的边界偏差。

3.2 先验与超参数 κ：无信息不等于无影响

很多人以为 α 全 1 就万事大吉，真正左右后验形状的是 κ。κ 很大时似然极度尖锐，后验几乎坍缩到 linprog 那个顶点附近，贝叶斯化白做；κ 很小时似然被压平，后验几乎等于先验，排序变成随机数。中间那段才是有效区间，判断办法是做敏感性分析：κ 取 0.1、1、10 各跑一遍，看 credal ranking 的结论是否稳定。稳定的那部分结论才是数据支持的，随 κ 翻脸的结论说明专家打分的信息量本来就不够。

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_BB=1、a_WW=1 是定义而不是观测，塞进似然等于给权重加了一个不该有的锚。第三，循环里两支的结构不对称：Best 那一支的分子是固定的 w_B，Worst 那一支的分母是固定的 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 = i+1:n if P(i,j) >= thr pairs{end+1} = sprintf('w%d > w%d (%.2f)', i, j, P(i,j)); elseif P(j,i) >= thr pairs{end+1} = 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 下出现了 A>B、B>C、C>A 的环，说明这三个准则在数据上根本分不开，该做的是回去检查打分或者合并准则，而不是换个随机种子重跑。

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_BB=1 和 a_WW=1 是定义值，不是专家给出的判断。把它们当数据放进似然，相当于额外加了两条“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

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