资讯动态

MATLAB贝叶斯决策模拟:从后验概率到最小风险分类

发布时间:2026/10/13 23:44:58 来源:尧图企业网站定制
1. 贝叶斯决策模拟从哪开始两类高斯数据与先验设置贝叶斯决策理论在模式识别、信号检测和机器学习里几乎是绕不开的基石但很多人在课堂上听得懂公式一上手用MATLAB模拟就卡住了。两年前我做某图像处理Demo时需要在两类数据之间做一个分类器原型验证当时想直接调现成库转头一想与其把分类器当黑盒调参不如自己把贝叶斯决策的过程完整模拟一遍既能看清后验概率的变化也能把最小错误率、最小风险的差别落到数值上对比。这个标题下要做的模拟核心就是把“给定观测x哪种类别最合理”这个问题变成可运行、可画图、可算错误率的实验。先明确模拟题目假设有两类样本分别服从二维高斯分布类别记作ω1和ω2。先验概率P(ω1)P(ω2)0.5类条件概率密度分别为p(x|ω1)和p(x|ω2)。用MATLAB生成训练数据、按贝叶斯公式计算后验概率、画决策边界、评估错误率。这是最基础的二元决策问题也是后续扩展到最小风险决策、三分类甚至连续观测的基础。为了让参数可调可复现我建议把均值向量和协方差矩阵作为脚本开头的显式变量% 类别1的均值与协方差 mu1 [1; 1]; S1 0.8 * eye(2); % 类别2的均值与协方差 mu2 [3; 3]; S2 0.8 * eye(2); % 先验 p_w1 0.5; p_w2 0.5;选均值差距为2、协方差为0.8I这个参数区间既能看出明显可分又不至于完全分开后验概率在重叠区域还有平滑过渡非常适合观察决策边界。如果均值差距拉到5以上两类基本不重叠模拟结果几乎和常识一致反而看不出贝叶斯决策的微妙之处。贝叶斯公式在这里的具体形式就是P(ωi|x) p(x|ωi) * P(ωi) / (p(x|ω1) * P(ω1) p(x|ω2) * P(ω2))其中p(x|ωi)用高斯密度函数计算。这里有一个初学者常忽略的数值经验对二维高斯密度取对数再计算后验比直接算指数更稳因为当x离均值较远时指数项会非常小直接算密度容易下溢。MATLAB里mvnpdf自带数值稳定处理但如果你想手写理解原理务必用对数域变换最后再还原。2. 网格化决策与决策边界可视化为什么不用测每个样本点模拟贝叶斯决策最容易想到的做法是在样本点上逐一计算后验概率、做出决策、统计错误率。这个流程本身没错但只做这一步会丢掉贝叶斯决策最有价值的信息——决策边界和风险曲面。我建议在样本点评估之外额外建立一个覆盖特征空间的网格在网格上计算后验概率再用等值线把边界画出来。这个网格化做法的好处在于它能显示完整分类面的形状而不只是看到几个点的分类结果。网格步长需要权衡步长太大边界锯齿严重步长太小计算量指数上升。以二维特征空间为例从-2到6的范围步长0.05就能得到160×160的网格对MATLAB来说是眨眼间的事[X1, X2] meshgrid(-2:0.05:6, -2:0.05:6); X_grid [X1(:) X2(:)]; % 把网格点排列成N×2矩阵计算网格上每个点的后验概率代码和样本点上是同一套函数。建议先封装一个函数避免在样本评估和网格评估之间复制粘贴function [post1, post2, decision] bayes_posterior(X, mu1, S1, mu2, S2, p_w1, p_w2) p_x_w1 mvnpdf(X, mu1, S1); p_x_w2 mvnpdf(X, mu2, S2); denom p_x_w1 * p_w1 p_x_w2 * p_w2; post1 p_x_w1 * p_w1 ./ denom; post2 p_x_w2 * p_w2 ./ denom; decision double(post1 post2); % 1表示判为类10表示判为类2 end关于mvnpdf的输入格式MATLAB要求均值是行向量、X是N×D矩阵、协方差是D×D矩阵。我最初犯过一次错把均值写成列向量[1;1]直接传进去结果报维度错误后来统一规定所有均值向量在函数内部转成行向量再使用这样外部传参再怎么折腾都不会出错。决策规则这里用的是后验概率比较也就是最小错误率决策当P(ω1|x) ≥ P(ω2|x)时判为ω1。看起来是在比较后验概率但代入贝叶斯公式会发现分母相同实际比较的是p(x|ω1)P(ω1)和p(x|ω2)P(ω2)。这就是为什么决策边界上的点满足p(x|ω1)P(ω1) p(x|ω2)P(ω2)对两个高斯分布取对数并化简决策边界其实是一个二次曲线两个协方差相等时退化为直线。这个推导值得自己手推一遍把高斯表达式的指数项展开令两者相等二次项抵消后得到关于x的线性方程。我在实验里取协方差相等所以边界是条直线这也是验证模拟正确性的一个天然标尺。网格化之后画图用contour画后验概率等值线用contour(f)或者contourc提取等值线坐标存下来还可以用fill在两条等值线之间填充颜色。我习惯的显示方案是散点图表示训练样本类1用圆圈、类2用叉号contour画出后验概率等于0.5的等值线作为决策边界再叠加一个背景热力图用imagesc或者pcolor显示后验概率的空间分布。这里值得强调一句决策边界是后验概率等于0.5的等值线而不是两类样本中间的垂直平分线。只有当先验相等且协方差矩阵相等时边界垂直平分两个均值点的连线一旦先验或者误判代价不对称边界会偏移。偏移的方向和幅度正是后面风险决策要重点讨论的内容。3. 从最小错误率到最小风险损失矩阵改写了决策规则光比较后验概率只回答了一个问题哪类的概率更大但在很多现实决策里两类错误造成的代价完全不一样。拿医学筛检打比方把健康人误判为患者后续检查成本高但后果相对可控把患者误判为健康人错过的可能是早期治疗窗口代价极大。这种情况下“哪类概率大就判哪类”就不再是合理策略需要引入损失矩阵。贝叶斯决策框架里定义损失λij为真实类别是ωj却判成ωi的代价。决策规则从“最大化后验概率”变成“最小化期望风险”——条件风险为R(αi|x) λi1 * P(ω1|x) λi2 * P(ω2|x)选择让条件风险最小的那个动作。等价地把待决类别的条件风险展开就会发现决策准则变成比较加权的似然比p(x|ω1) / p(x|ω2) 与 [λ12 * P(ω2)] / [λ21 * P(ω1)] 比较这里λ12表示把类1误判为类2的代价λ21表示把类2误判为类1的代价。这个比值关系非常重要它告诉你误判代价的提升等价于推着决策边界向另一边移动让决策者更倾向于避开代价高的错误。在MATLAB里实现这个规则不需要另起炉灶。沿用前面的后验概率函数加一个损失矩阵参数然后比较条件风险% 损失矩阵L(i,j) 表示真实类别为j判为i的损失 L [0 1; 1.5 0]; risk_alpha1 L(1,1) * post1 L(1,2) * post2; % 判为类1的条件风险 risk_alpha2 L(2,1) * post1 L(2,2) * post2; % 判为类2的条件风险 decision_risk double(risk_alpha1 risk_alpha2);我用L(2,1)设置为1.5、L(1,2)保持1意味着把类1的样本误判为类2比把类2误判为类1的代价低一半。直觉上决策器应该更“愿意”把样本判给类1边界应向类2方向移动。实验也确实符合决策边界从原来的x1x24这条垂直线整体朝类2方向偏移了大约0.3个单位。为什么偏移不是0.5个单位因为损失比和先验比是联合作用的。我们来算一笔细账P(ω1)P(ω2)0.5λ21/λ121.5阈值就是λ12P(ω2) / (λ21P(ω1)) 1/1.5 ≈ 0.667。决策边界发生在似然比等于0.667处而不是1。在两个等协方差高斯的情形下对数似然比是x的线性函数阈值从0变成ln(0.667)-0.405等价于边界平行移动距离 0.405 / ||w||这里的w是线性判别方向的法向量。模长||w||大约1.25所以移动约0.32个单位和实验观察吻合。这个仿真过程强烈建议读者亲手跑一遍。你会在图上看到一个非常直观的现象先验和损失的不对称最终都归结为决策阈值的平移而边界形状在等协方差时是直线基本不变。这个洞察对理解许多分类器里的偏置项bias很有帮助——很多线性分类器的偏置本质上就是在调整这个阈值只是有人理解成了“模型参数”有人理解成了“贝叶斯阈值”。4. 蒙特卡洛评估与贝叶斯风险错误率数值怎么算才算对建造好了决策器接下来要面对一个更麻烦的问题这个决策器到底好不好最直接的指标是错误率而错误率有两种算法对应两种视角。第一种是理论视角。给定决策边界把特征空间分成R1和R2两个区域错误率是P(error) ∫_{R2} p(x|ω1)P(ω1)dx ∫_{R1} p(x|ω2)P(ω2)dx大多数时候这个积分没有解析解但对等协方差高斯、线性边界的特例可以用一维正态的累积分布函数表示。MATLAB里用normcdf直接算% 在投影方向w上类条件分布是一维高斯 % 阈值 t 由决策边界在w方向上的截距决定 error_theory P_w1 * normcdf((t - mu_proj1)/sigma_proj) ... P_w2 * (1 - normcdf((t - mu_proj2)/sigma_proj));第二种是实验视角也就是蒙特卡洛模拟。按先验比例生成大量测试样本比如每类各5000个用决策器判定统计误判比例。两种视角互为印证如果理论值和模拟值差得远说明边界画错了或者数据生成和预期不一致这个对照本身就是调试工具。我在实践中的经验是两者应吻合到小数点后两位一般在0.01以内都算正常。蒙特卡洛统计错误率的代码骨架rng(42); % 固定随机种子保证可复现 N_test 5000; % 按先验比例生成样本 n1 round(N_test * p_w1); n2 N_test - n1; X_test1 mvnrnd(mu1, S1, n1); X_test2 mvnrnd(mu2, S2, n2); [~, ~, dec1] bayes_posterior(X_test1, mu1, S1, mu2, S2, p_w1, p_w2); [~, ~, dec2] bayes_posterior(X_test2, mu1, S1, mu2, S2, p_w1, p_w2); err1 sum(dec1 0) / n1; % 类1被误判为类2的比例 err2 sum(dec2 1) / n2; % 类2被误判为类1的比例 error_total err1 * p_w1 err2 * p_w2;这里的rng(42)值得多说一句。很多初学者做蒙特卡洛时忽略了随机种子每次运行结果都变写报告时不知道拿哪组数字。固定种子后只要参数不变实验结果完全可复现这在做对比实验时尤其重要。贝叶斯风险比错误率更进一步错误率隐含假设“两种错误的代价都是1”而贝叶斯风险把损失矩阵的代价乘进去R(δ) λ12 * P(w1) * err1 λ21 * P(w2) * err2在最小风险实验里用带损失矩阵的决策规则计算出来的贝叶斯风险一定不高于用最小错误率规则直接算的风险。这个结论在理论上必然成立但模拟中如果不加验证容易因为代码错误得到荒谬的相反结果。我做了一个快速验证用最小错误率规则跑一遍风险约为0.92因为错误率虽然较低但误判类1样本的损失被1.5倍放大换成最小风险规则后风险降到0.83左右。数字本身不重要重要的是在同一组测试集上两条决策规则的风险值能分出高低逻辑闭环。由此引出一个实操层面的注意点评估决策器的指标必须和优化目标一致。如果你用最小风险规则训练决策器却只用错误率评估很可能会得出“最小风险规则不如最小错误率规则”的片面结论——因为前者在错误率的定义下确实略高但它在真正关心的代价维度上更优。模拟实验里能同时追踪多个指标时不要只看一个。5. 可视化呈现的进阶技巧把后验概率和边界画成有说服力的图表模拟的最后一步是输出结果而MATLAB可视化恰恰是很多人的短板。只画一堆散点和一条边界线报告里看起来单薄如果能把后验概率的分布、样本的落区、边界的形态一起呈现出来效果完全不同。我的标准画法分四个子图第一个子图散布样本点用不同的标记和颜色区分真实类别gscatter可以一句搞定gscatter(X(:,1), X(:,2), Y, rb, xo);第二个子图用contourf绘制后验概率P(ω1|x)的填色图。颜色从深到浅代表后验从高到低边界正好在颜色过渡的中央。这里我推荐一个技巧把colormap设置为蓝色到红色渐变这样视觉上天然产生“一类蓝区、一类红区”的观感边界清晰可见contourf(X1, X2, reshape(post1, size(X1)), [0.1 0.3 0.5 0.7 0.9]); colormap([linspace(0.2,1,64) ones(64,1)*0.2... ]); % 具体用bluewhitered更省事第三个子图叠加决策边界和测试样本的分类错误点。错误点用黑色方块标出来会让错误率的来源一目了然比只给一个总数更有说服力mis_idx (Y ~ pred); plot(X_test(mis_idx,1), X_test(mis_idx,2), ks, MarkerFaceColor, k);第四个子图画出理论边界和实际边界的对比。如果两个边界重合说明模拟实现正确不重合就倒回去找代码问题。这个子图是我个人最喜欢的调试工具因为它的反馈信号非常明确——重合就是全对不重合就是有问题不需要猜。画图时还有一个小坑contourf的输入矩阵必须是二维的而我们在网格上计算的后验结果是拉直的长向量需要先reshape回网格形状。很多初学者在这一步报错就是因为忘了reshape。我的经验是在函数返回后立即reshape并存入同名变量比如post1_grid reshape(post1, size(X1))后面所有画图都用这个网格版本避免混淆。可视化里另一个值得做的实验是“后验概率随样本位置的连续变化”在特征空间里取一条穿越边界的线段把后验概率沿这条线画成曲线。你会发现远离边界时后验趋近于1或0靠近边界时后验平滑过渡曲线斜率由协方差矩阵决定。协方差越小曲线越陡决策越自信协方差越大曲线越缓边界附近越“含糊”。这种图对理解置信度的空间分布非常有帮助也是把贝叶斯决策和“分类器置信度校准”联系起来的一个桥梁。6. 扩展实验敏感损失下的决策偏移与三分类问题前面的实验都是两类、二维、等协方差高斯条件是刻意简化的。实际应用中很少这么规整所以至少再跑两组扩展检验模拟框架的普适性。第一组扩展验证损失矩阵对决策边界的平移效应。取λ21从1.0扫到3.0步长0.1对每个值重新计算决策边界和贝叶斯风险。把边界偏移量画成随λ21变化的曲线会看到一条近似线性的递增关系。我实测过λ21从1到3边界偏移从0增加到0.8左右而错误率变化却很小因为总错误率本来就在0.1左右徘徊加上样本重叠区域有限边际影响被稀释。这组实验说明了一个很容易被忽视的事实损失矩阵变化主要影响风险的构成而错误率并不是最敏感的指标。第二组扩展把两类推广到三类。三分类的贝叶斯决策不再是“找一个阈值”而是划分三个区域每个区域对应后验概率最大的类。决策边界不再是一条直线而是可能交汇于一点的组合边界甚至可能是曲线。MATLAB代码改写的关键在于把后验概率函数推广成循环结构function post bayes_posterior_k(X, mus, Sigmas, priors) K length(priors); n size(X, 1); post zeros(n, K); for k 1:K post(:, k) mvnpdf(X, mus{k}, Sigmas{k}) * priors(k); end post post ./ sum(post, 2); end决策规则变成[max_val, decision] max(post, [], 2)边界用contourf画每个类别的获胜区域。三分类实验里一个有意思的现象是即使三类先验完全对称边界也不一定是等距的——因为协方差矩阵不同边界会向高方差类别方向凹陷这正是贝叶斯决策对“不确定性更大”区域给予更多宽容的体现。扩展实验做下来我个人最大的体会是贝叶斯决策的框架并不复杂但每增加一个变量代价不对称、类别数增加、协方差不等决策结果都会发生微妙而可解释的变化。如果只记公式不模拟这些细微变化很难真正在直觉中扎根。而MATLAB的矩阵运算把实现成本压得很低配合可视化几乎是最适合做这类实验的环境。7. 从模拟回归应用场景贝叶斯决策和匹配滤波的关系贝叶斯决策模拟做完以后如果你只把它当成一个分类作业题就浪费了这套思路的潜力。它在信号处理里的经典对应是匹配滤波而这个映射在模拟里可以顺手验证。考虑一个简单场景发送端以等概率发送两个确定信号s1(t)和s2(t)接收信号x(t)受到高斯白噪声污染。目标是判断发送的是哪个信号。把连续波形在观测区间内离散成N个点观测向量x就是N维随机向量两类假设分别是x ~ N(s1, σ²I)和x ~ N(s2, σ²I)。这和我们前文的二维高斯分类在数学上完全同构只是维度从2升到N。在N维情形下贝叶斯决策的充分统计量是线性相关运算即把接收信号和[s1-s2]做内积再与阈值比较。这就是匹配滤波的本质匹配滤波是高斯噪声下贝叶斯决策的直接推论。用MATLAB模拟时不必真的穷举所有N维网格只需要在测试样本上计算判决统计量并统计错误率理论错误率用Q函数或者normcdf算SNR sum((s1 - s2).^2) / sigma^2; P_error qfunc(sqrt(SNR/4)); % 等先验、等代价时的最小错误率我在模拟里常用s11的方波和s2-1的方波噪声方差从0.1逐步增加到10画出错误率随SNR变化的曲线理论和模拟两条曲线会几乎完全重合。这个实验对理解“为什么相干接收机能达到最优性能”特别有帮助——不是因为它用了什么神秘算法只是因为它实现了贝叶斯最优判决。如果读者有通信系统仿真背景可以把这段当成贝叶斯决策和信号检测之间的桥梁。很多讲模式识别的教材不提匹配滤波讲通信的教材又不推导贝叶斯决策但数学上它们是同一件事。自己做一次模拟把它们连起来后面再遇到相关问题时会有一种“原来如此”的通透感。8. 踩坑记录贝叶斯模拟中五个容易翻车的地方模拟过程中踩过的坑整理出来供大家参考。第一个坑是协方差矩阵非正定。很多人喜欢手动输入数据生成协方差比如S[0.8 0.6; 0.7 0.9]看起来差不多但对称性不满足mvnpdf直接报错或者结果荒谬。解决方法是构造时用S[a b; b c]的形式保证对称再用一个对称正定矩阵的构造技巧A randn(2,2); S AA 0.01eye(2)。后者保证正定适合快速实验。第二个坑是数据生成顺序和先验不一致。比如先验设置为P(ω1)0.7但生成测试样本时两类各生成5000个那么蒙特卡洛评估的总错误率公式里乘的权重就和实际样本比例对不上计算结果会失真。正确的做法是严格按照先验比例生成样本或者用样本比例作为权重。我建议固定总样本数然后按先验比例取整。第三个坑是网格范围不够大。当协方差矩阵特征值较大时高斯分布延展很广如果网格范围只覆盖均值附近后验概率在边缘处可能仍然高于0.5因为另一类更远决策边界被截断。我的经验是先看两类均值的位置向外扩3到4倍标准差再看边界是否完全包含在网格内。第四个坑是决策边界的提取方式。用contour直接画等值线简单但想得到边界对应的精确解析位置建议改用解析式计算。等协方差时边界是线性方程(x - x0)^T * S^{-1} * (mu1 - mu2) 0其中x0 (mu1mu2)/2。用这个公式画出的边界直线和contour的0.5等值线应该完全重合可以作为正确性的交叉验证。第五个坑是后验概率计算时忽略了先验。有些代码写出来是p_x_w1和p_x_w2直接比较没有乘先验这在先验相等时没问题但一旦先验变化决策边界完全错位。这类bug很隐蔽因为结果“看起来也合理”。我通常会设计一个不对称先验的测试用例专门检查先验设为P(w1)0.9、P(w2)0.1看边界是否向类2方向移动如果没有代码一定漏了先验乘法。下面是快速诊断对照表现象最可能原因检查方式理论错误率和模拟错误率不符边界解析式或蒙特卡洛权重有误打印边界方程系数和阈值逐步核对决策边界完全偏到一侧先验或损失矩阵没有参与计算改为对称参数确认边界回到几何中线网格上后验概率全为0或1网格范围太小或协方差接近奇异扩大网格检查协方差矩阵条件数运行时报错维度不一致均值向量行/列方向错误统一在函数内转成行向量每次运行结果差异大未固定随机种子脚本开头加rng(固定值)9. 参数速查表与完整模拟流程总结为了方便直接照着跑我把模拟中用到的关键参数和函数整理成下面的速查表。项目推荐取值/写法说明均值向量mu1[1;1], mu2[3;3]可换但不要离太近否则几乎不可分协方差矩阵S10.8eye(2), S20.8eye(2)等协方差时边界是直线便于验证先验p_w10.5, p_w20.5换成0.7/0.3观察边界平移损失矩阵L[0 1; 1.5 0]第二行第一列大于1判决偏向类1网格范围-2:0.05:6覆盖两类均值外扩至少2个标准差测试样本量每类5000足够小和理论对比稳定性好随机种子rng(42)固定后可复现后验概率函数前文给出的bayes_posterior建议封装复用理论错误率normcdf 或 qfunc仅等协方差线性边界时有解析式整个模拟流程遵循一个清晰的闭环设计问题、定义参数、生成数据、计算后验、画决策边界、蒙卡评估、对照理论、调整参数再验证。每走完一个循环对这个决策问题的理解就深一层。这一套框架不只适用于高斯分布假设。把mvnpdf替换成其他概率密度的评估函数比如泊松、指数或者自定义核密度整个决策框架照常运转因为贝叶斯决策的核心不依赖分布形式只依赖后验概率的计算。这是这套模拟方法最值得保留的通用性。最后说一个我用了很久的小习惯每次模拟都在脚本开头写一段注释记录这次实验要回答什么问题、关键参数是怎么选出来的。两个星期后回看时这段注释比任何代码都宝贵因为参数背后的决策动机已经融进了注释里。模拟贝叶斯决策表面上是在做数值实验本质上是在训练一种“用代价和概率思考问题”的习惯。这个习惯才是比任何一行MATLAB代码都值钱的东西。

读完文章,也想定制专属网站?

尧图设计师 24 小时内与您沟通定制方案

免费获取报价 →
↑