简介这是一套基于MATLAB实现的多模型集成BMA方法资源面向气象、环境、水文等领域的科研人员与中高级学习者旨在解决单一模型预测不确定性高、精度不稳定等问题。资源以贝叶斯模型平均为理论内核在备选模型集合中通过验证集表现计算各模型权重再按权重执行加权平均从而降低整体预测风险、提升稳健性该策略对样本有限、模型差异明显的场景尤其有效。压缩包共13个文件其中8个M脚本覆盖数据预处理、期望最大化求解、BMA权重计算、模型拟合、散点与密度图绘制及误差对比等完整流程2个MAT数据集提供预测与观测样例便于直接运行验证另有2个TXT说明文档和1个ZIP扩展包整体约5.23MB模块划分清晰从预测数据出发经核心脚本分配权重后可输出融合预测并与观测数据比对。目前已有842人学习下载适合希望快速上手多模型集成方法并借助完整代码理解算法细节的读者。借助这些脚本可系统掌握BMA的权重推导与程序实现复用可视化工具检验模型性能为后续气候或环境预测研究搭建可扩展的实验基底。1. 多模型集成的老问题BMA 不是从模型里挑一个而是给每家的预测配权重气候预测里最常遇到的一个尴尬局面是你手里有好几个模型的结果可它们互相矛盾。有的系统性偏低有的在高值段发散还有的在某些年份表现特别好、另一些年份一塌糊涂。强行选一个最好的等于把其他模型里有效的信息全部扔掉直接做简单平均又把每个模型的偏差也一起平均进去了。多模型集成 BMA 解决的就是这个怎么分配信任的问题先对每个模型的输出做偏差校正再通过似然函数在训练期上计算每个模型的后验权重最后对校准后的预测做加权平均。它适合所有需要把多套预测结果整合成一个稳定预测的场景尤其是气候、水文和环境模拟这类模型各有适用区间、没人知道哪个绝对正确的现实局面。这套资源是用 MATLAB 实现的多模型集成 BMA 完整流程脚本结构清晰你拿去改数据路径就能跑。2. BMA 数学模型条件分布、EM 迭代、偏差校正三块决定权重质量在动手碰代码之前先花一页纸把 BMA 的内部逻辑说清楚。这一章的三个小节对应项目里的三个核心脚本fit.m 做偏差校正EM_new.m 做迭代求权重BMA_new.m 把它们串起来。理解它们的先后顺序比你直接跑通脚本重要得多。2.1 核心公式与两个约束预测分布是权重的混合分布假设有 K 个模型每个模型在第 n 个时间点给出一个预测值 f_{n,k}观测值是 y_n。BMA 的预测分布是每个模型条件分布的加权混合p(y_n | f_{n,1}, ..., f_{n,K}) sum_k w_k * g_k(y_n | f_{n,k})其中 g_k 是模型 k 在校正之后的条件密度函数最常见的选择是正态分布。权重 w_k 有两个约束w_k ≥ 0且 sum_k w_k 1。这两个约束保证了混合结果是一个合法的概率分布也是 BMA 与一般回归加权最大的区别——权重永远非负不会出现某个模型权重是 -0.2这种没法解释的结果。这里有一个新手常有的疑问为什么不直接用每个模型的训练期 RMSE 换算成权重RMSE 只能反映确定性预测的误差但 BMA 关心的是给定这个模型时观测落在哪的概率最大它同时考虑了偏差和不确定性。两个模型 RMSE 相同一个残差是正态分布一个残差重尾BMA 会倾向于前者。这一点在后续 EM 的实现里会体现在正态分布的对数似然计算上。具体到代码层面g_k 的参数一般表达成 N(a_k b_k * f_{n,k}, sigma_k^2)。a_k 和 b_k 来自 fit.m 的线性回归拟合sigma_k 是残差标准差。每个模型在参与权重计算之前都会先被单独校准一遍这是 BMA 能应对气候模式系统性偏差的关键设计。2.2 EM 算法E 步算责任M 步更新权重直接最大化 BMA 的对数似然很困难因为隐变量——每个时刻到底是哪个模型真正产生了观测——是不可观测的。标准解法是 EM 算法它把最优化问题拆成交替两步。E 步对每个模型 k 和每个时间点 n计算后验责任 gamma_{n,k}表示第 n 个时刻模型 k 承担解释观测的责任gamma_{n,k} w_k * p(y_n | f_{n,k}) / sum_j w_j * p(y_n | f_{n,j})分母就是当前权重下的混合密度分子是模型 k 的加权贡献。计算时要注意数值稳定性对数密度直接 exp 很容易溢出常见做法是先把对数密度减去该行的最大值再 exp避免中间结果出现 inf 或 NaN。M 步更新权重为责任的时间平均w_k^(new) (1/N) * sum_n gamma_{n,k}这个更新就是 EM 的极大化步骤表达式简单到不敢信你只需要把 E 步算出来的责任矩阵按列求均值。E 步和 M 步交替进行直到对数似然的变化小于某个阈值我一般取 1e-4或达到最大迭代次数200 次足够。具体到正态假设下对数密度写出来是log p(y_n | f_{n,k}) -0.5 * ((y_n - mu_{n,k}) / sigma_k)^2 - log(sigma_k) - 0.5 * log(2*pi)其中 mu_{n,k} a_k b_k * f_{n,k}。sigma_k 出现在分母上这保证了残差更小的模型在对数似然上天然占优但注意这种占优是相对概率层面的而不是直接把 sigma 的倒数当作权重。一个常被忽略的细节EM 对初始权重并不免疫。如果初始权重把某个模型压到接近 0后续迭代很难摆脱这个困境。常见做法是从均匀权重 w_k 1/K 出发或者对每个模型单独做一个以它为主的初始化最后比较哪个收敛点的对数似然最高。这一点在避坑章里还会展开。2.3 偏差校正先行为什么 BMA 不能拿原始预测直接算权重很多第一次接触 BMA 的人会跳过偏差校正直接把模型原始输出代入正态密度函数。这在大多数气候数据上是错误的。原因很直接气候模式输出往往带明显的系统性偏移整体偏高一度或偏低几毫米都很常见原始残差方差巨大BMA 算出的权重会被方差差异主导而不是被模型的可校正性主导。说白了一个偏差大但趋势灵敏的模型在校准之后可能是最有价值的但如果直接用原始数据算它会被误判成最差的。fit.m 在项目里承担的就是这一步。它的输入是两个长度相同的向量观测 y 和某个模型的预测 f_k输出是线性回归的系数 a_k、b_k 以及残差标准差 sigma_k。拟合完成之后参与 BMA 的预测值变为 a_k b_k * f_k残差计算公式变为 y_n - (a_k b_k * f_{n,k})。这一步还有一个容易被忽略的细节偏差校正的回归应该在训练集上完成然后用同一组 a_k、b_k 去转换验证集的预测值千万不能把验证集也放进回归拟合。否则验证阶段的似然是被污染的权重会虚高。这个操作和交叉验证的先拟合、再评估原则完全一致。在实际项目里偏差校正的有无对权重结果的影响是肉眼可见的。某次我用气候模式数据做春季降水集成跳过 fit.m 直接跑 EM结果权重几乎全压在方差最小的那个模型上加上偏差校正之后权重分布才变得符合直觉——表现稳定且偏差小的模型拿到大头一个有偏差但趋势敏感的模式也拿到了接近 0.2 的份额。所以不管你的数据是从哪来的这条流程都建议不要省。3. 代码包地图与数据准备先分清核心脚本和工具函数再对齐观测与预测这个压缩包里文件不少但实际上分工很清晰。如果一上来就双击 BMA_new.m 跑很可能因为数据格式不符直接在 load 阶段报错。我建议先花几分钟把文件分成几类心里有个地图再动手。3.1 文件清单与分工四条逻辑主线整个包按功能可以分成四组。第一组是核心算法链。BMA_new.m 是主入口负责读数据、调用 EM、计算权重、最后输出集成预测EM_new.m 是 EM 迭代的函数体一般接收观测值、每个模型的校准预测序列、初始权重和迭代参数返回最终权重和每个时刻的后验责任矩阵fit.m 做偏差校正线性回归输出 a、b、sigma 三个参数。这一组是你需要真正读懂的全部。第二组是数据文件。CN_predata.mat、CN_obsdata.mat、data.txt。前两个是 mat 格式的观测与预测序列CN 前缀大概率表示某个区域尺度的数据集predata 是各模型的预测序列obsdata 是对应的观测序列。data.txt 是同一份数据的文本版方便你不启动 MATLAB 也能先看看数据大概长什么样。第三组是可视化模块。scatplot.m、dataDensity.m、DataDensityPlot.m、TestDataDensity.m。dataDensity.m 负责计算样本点的密度估计scatplot.m 把密度渲染成散点颜色DataDensityPlot.m 是封装好的绘图脚本TestDataDensity.m 是可视化模块的测试入口。它的作用不是参与 BMA 计算而是帮你评估集成预测的质量。第四组是扩展和说明。BMA_IPCC_yr.m 是把 BMA 按年度循环处理多模式集成的版本适合逐年份输出权重license.txt 是许可证DataDensity.zip 大概率是密度绘图相关的补充文件。一个常见误区是把 DataDensityPlot.m 当成 BMA 主程序。从文件名和依赖关系看它只调用 dataDensity 和 scatplot不调用 EM_new.m所以它只是数据展示层权重计算完全与它无关。3.2 mat 数据与 txt 数据的结构先看列数再看年限拿到数据的第一步不是跑代码是确认维度。在 MATLAB 命令窗口执行load(CN_predata.mat); load(CN_obsdata.mat); whoswhos 会列出工作区里每个变量的名字、尺寸和类型。常见的结构是 predata 是一个 T×K 矩阵T 是时间长度K 是模型数量obsdata 是 T×1 列向量。但也存在另一种常见布局数据被存成结构体比如 S_pred.years、S_pred.models这时候需要先看字段名再取数。接下来对齐时间。气候数据最怕的就是年份错位观测从早期开始某个模型可能晚几年才开始。建议先检查每列是否都有完整覆盖obs CN_obsdata; % 假设是 T×1 pred CN_predata; % 假设是 T×K size(obs) size(pred)如果发现 pred 行数比 obs 多说明时间范围不一致需要按年份截取。data.txt 如果存在通常每一行是一条记录前几列是年份、月份后面是各模型预测和观测可以用 readtable 快速读进来核对。3.3 数据对齐与清洗清掉 NaN 再谈权重在把数据交给 EM 之前有两件必须做的事时间对齐和非有限值清洗。气候观测里有缺测是很正常的而缺测会直接让 EM 的对数似然变成 NaN然后整片权重跟着失效。% 用 isfinite 找出所有列都是有限值的行 good isfinite(pred) isfinite(obs); keep all(good, 2); % 行方向求与得到 T×1 逻辑向量 pred_clean pred(keep, :); obs_clean obs(keep); % 如果数据带年份向量先确认两端年份一致 years_obs obs_year; % 观测年份 years_pred pred_year; % 预测年份假设和观测相同需要注意 isfinite 是逐元素的all(..., 2) 是按行判断这一行所有模型都没有缺测。这句话的意思是只要某个模型在某一年缺测整行都会从数据集里剔除。这样处理的代价是损失一些样本但保证了似然函数的完整性从稳定性角度更值。如果缺测的模型只有一两个更精细的做法是把缺测行的模型列用该模型的多年平均值插补而不是整行删除这个在样本量本来就小的站点上尤其推荐。我个人的习惯是先看缺测比例低于 5% 用删除高于 5% 用均值插补再高就需要考虑是不是数据下载环节出了问题。对齐和清洗做完主流程才值得继续。这里顺带提醒清洗后的行数要记录下来最后画图和算 RMSE 时也要用同一组索引否则集成预测和观测对不齐可视化的散点图会出现错位。4. 跑通 BMA_new.m从载入数据到输出权重的完整代码路径这一章是完整抄作业的部分。以下代码的结构按照常用的 BMA 气象集成实现方式补全与原包内脚本的功能一一对应。如果你的版本在细节上有差异重点看函数之间的数据流转主脚本把数据传给 fit.mfit.m 返回校准参数EM_new.m 用校准后的序列迭代权重最后主脚本做加权集成。4.1 载入与初始化先跑 fit.m把每个模型的校准参数存下来%% 载入数据 load(CN_predata.mat); load(CN_obsdata.mat); % 假设 pred 是 T×K 矩阵obs 是 T×1 列向量 T size(pred, 1); K size(pred, 2); % 清洗保留所有模型和观测都非 NaN 的行 keep all(isfinite(pred), 2) isfinite(obs); pred pred(keep, :); obs obs(keep); N sum(keep); %% 偏差校正对每个模型做线性拟合 y ~ a b*x mu zeros(N, K); % 校准后的预测均值 sigma zeros(K, 1); % 每个模型的残差标准差 for k 1:K [a_k, b_k, s_k] fit_data(obs, pred(:, k)); % 对应包里的 fit.m mu(:, k) a_k b_k * pred(:, k); sigma(k) s_k; end % 权重初值均匀分布 w ones(K, 1) / K;fit_data 内部一般就是用 polyfit 或 regress 拟合一次线性方程然后把残差的标准差取出来。注意 sigma 的估计也可以用自由度修正版本即除以 N-2 而不是 N小样本时误差更小这个细节在 fit.m 里通常已经处理。这里还有一个值得关注的参数如果某个模型预测的范围和观测不在同一量级比如观测是 0~100 mm 降水某个模型输出却是 0~1 的标准化值fit.m 的回归系数会自动吸收这个尺度差异但数值稳定性会变差。我一般在跑之前会对预测量做均值方差标准化输出权重之后再还原回原始尺度这样 EM 的迭代过程不会因为数值尺度问题震荡。4.2 EM 主循环对数似然监控与收敛判断%% EM 迭代E 步与 M 步交替 max_iter 200; tol 1e-4; loglik_old -inf; for it 1:max_iter % E 步计算每个模型在每个时刻的对数密度 % mu 是 N×Ksigma 是 K×1obs 是 N×1 log_dens -0.5 * ((obs - mu) ./ sigma) .^ 2 ... - log(sigma) - 0.5 * log(2 * pi); % 数值稳定性每行减去该行最大值避免 exp 溢出 max_log max(log_dens, [], 2); log_dens_shift log_dens - max_log; % 未归一化的责任分子部分 weighted w .* exp(log_dens_shift); % 归一化den 是每行混合密度gamma 是后验责任 den sum(weighted, 2); gamma weighted ./ den; % N×K 后验责任矩阵 % M 步权重 责任矩阵的列均值 w_new mean(gamma, 1); % 计算当前对数似然 loglik sum(max_log log(den)); % 收敛判断 if it 1 abs((loglik - loglik_old) / abs(loglik_old)) tol fprintf(EM 收敛于第 %d 次迭代loglik%.6f\n, it, loglik); break; end loglik_old loglik; w w_new; end w w / sum(w); % 确保权重归一化几个关键点拆开说。第一log_dens 的表达式就是正态分布的对数密度mu 和 sigma 已经在前一步被 fit.m 校准好了所以这里不需要再动偏差矫正。第二max_log 的减法是为了防止 exp 溢出这是 EM 实现里的常规操作去掉这行你的脚本可能在某个模型残差方差异常小时直接算出 inf。第三收敛条件用的是相对变化量而不是绝对变化量因为对数似然本身的量级会随数据量增大而变大绝对阈值不适合所有数据集。这里我不建议把收敛条件设得太松。1e-4 是经验值如果发现权重在连续迭代间反复横跳先检查是不是有两个模型预测高度相关、责任矩阵几乎相同这种情况下需要额外加一个很小的正则项比如权重加上 1e-6 再归一化帮助迭代稳定下来。4.3 输出与集成预测从权重到最终结果得到权重 w 后BMA 的确定性预测就是每个模型校准预测的加权平均%% 集成预测 y_bma mu * w; % N×1 加权平均预测量 % 概率预测蒙特卡洛抽样生成预测分布 M 5000; % 抽样次数 sim zeros(N, M); for k 1:K nk round(w(k) * M); % 该模型应分的样本数 if nk 0 start_col sum(round(w(1:k-1) * M)) 1; end_col sum(round(w(1:k) * M)); sim(:, start_col:end_col) mu(:, k) randn(N, nk) * sigma(k); end end这一段逻辑上更直观的做法是每个时刻以权重 w 为概率从 K 个模型中挑一个然后从该模型的分布 N(mu, sigma^2) 中抽一个样本。上面的代码为了保证每个模型被抽到的次数严格等于 w*M使用了近似分层抽样。实际使用时用 randsample 或者简单的离散抽样都行结果差异不大。最后输出的 y_bma 是你真正用来画图和评估的集成预测。注意它用的是校准后的 mu 而不是原始模型输出所以在和观测画散点图之前记得确认 y_bma 和 obs 的长度都是 N。做完这一步可以顺手把 w、y_bma 和 keep 索引保存成 mat 文件后面做交叉验证和敏感性分析时直接复用。5. 避坑专题BMA 权重计算中最常见的五个翻车点这一章是实打实的教训合集。每一条我都踩过或者看同行踩过按现象、原因、解决写清楚方便你对照排查。5.1 权重结果对初始值敏感EM 收敛到局部最优现象只改了 w 的初值迭代出来的权重分布差得很远有时模型 A 占 0.6有时模型 B 占 0.5两组结果的对数似然都很高你不知道该信哪个。原因EM 算法只能保证收敛到局部最优BMA 的对数似然面在很多模型强相关时会出现不止一个局部峰。解决跑多组初始值。最稳的做法是三个起点均匀权重、每个模型单独为主w_k0.9其余均分、随机权重分别跑完 EM 后比较对数似然取最大的一组作为最终结果。这一步计算量不大K 个模型跑 3~5 组初始值也就几秒到几十秒。5.2 预测数据里有 NaN权重整片变 NaN现象跑完 EM 后 w 全是 NaN或者某个模型的权重直接消失后面的集成预测 y_bma 也是 NaN。原因某个模型的回归残差里有缺测导致 sigma(k) 是 NaN进而在 log_dens 计算时整列变成 NaN责任矩阵随之全 NaN。解决进入 EM 前强制做 isfinite 过滤并且把过滤后的行数打印出来确认每个模型实际参与计算的有效样本量。如果缺测是零散几条删除整行可以接受如果某一列缺测超过 20%建议检查数据源而不是硬着头皮跑。5.3 交叉验证切分切错方向权重被高估现象用随机划分的 70% 训练、30% 验证BMA 在验证集上的 RMSE 比简单平均好一大截但实际业务里复现不出同样的效果。原因气候和水文数据是强自相关时间序列相邻年份的观测高度关联随机抽样会把几乎重复的年份同时放进训练和验证等于让模型变相看到了答案。解决时间序列必须用时间块交叉验证。常见做法是前 20 年训练、后 10 年验证或者滚动窗口保证训练期完全在验证期之前。你可以用循环把数据按年份切块逐块输出权重最后看权重是否在块间漂移。5.4 权重集中在单一模型别的模型权重接近 0现象训练期稍长一点权重就全部压到 RMSE 最小的模型上其他模型权重只有 0.01 甚至更小BMA 退化成选模型而不是集成模型。原因样本量较大时似然函数对残差方差的差距非常敏感一个略好的模型会吃掉所有概率质量。解决回到 BMA 的贝叶斯本质给权重加一个 Dirichlet 先验或伪计数。具体做法是在 M 步更新时改成 w_k (sum(gamma(:,k)) alpha) / (N K*alpha)alpha 取 0.01~0.1 就能明显缓解权重两极分化又不会过度改变结果。5.5 跳过偏差校正权重给到平均最接近而不是校准最准的模型现象直接用原始模型输出跑 EM算出来的权重反直觉——某个显然有系统偏差的模型拿到了很高权重。原因系统性偏差没有消除前残差方差会偏大BMA 的似然把所有模型的差异都算成了不确定性权重因此失真。解决永远先走 fit.m 的线性校准流程。如果不确定校准是否有效可以输出每个模型的残差标准差对比一下校准后的 sigma 通常显著小于原始残差BMA 的权重才会落在真正可信的模型上。6. 进阶验证BMA_IPCC_yr.m 多场景跑法与密度图诊断6.1 多场景跑法BMA_IPCC_yr.m 在包里的定位是按年循环的多模式集成版本适合把数据拆成逐年窗口跑 BMA然后对比各年权重。跑法上我通常是先像第 4 章那样完成一次全时段的权重计算再单独跑三个子时段前期、中期、近期。如果权重在三个子时段中明显漂移说明模型的可信度随时间变化这时直接做全时段集成反而不合适应该把时段切分为两个独立窗口分别出预测。6.2 密度图诊断DataDensityPlot.m 配合 scatplot.m 怎么看打开 TestDataDensity.m 跑一遍它读取一组预测-观测对在散点图上叠加密度色。读图时重点看三个位置主对角线附近颜色是否集中、离线区是否有明显的高密度块、远离对角线的离散点数量。如果 y_bma 和 obs 的散点密集贴合对角线说明集成预测的校准效果良好反之如果出现平行于对角线的偏移带就要回去检查 fit.m 的 a_k 是否被错误地设成了 0。6.3 与简单平均的对比评估 BMA 是否真的比简单平均好不能只看 RMSE。通常至少看两个指标确定性预测的 RMSE 和概率预测的 CRPS。经验上 BMA 的 RMSE 与简单平均相差不大但 CRPS 会明显更优因为 BMA 输出的是完整概率分布而不只是一个点。指标简单平均BMA 加权集成RMSE通常是基线略优或持平CRPS无概率信息公平性差通常明显更优权重可解释性无可解释、可检验漂移这套资源里 BMA_IPCC_yr.m 和密度图脚本合起来已经够做一个完整的多模型集成评估流程。你自己跑的时候建议把多初始值的对数似然曲线、三个子时段的权重变化、密度图这三样东西放在同一个报告页里互相印证比单看一个数字可靠得多。我印象最深的是某个降水集成项目里单纯跑默认参数得到的权重分布还算合理一换成多初始值扫描才发现默认起点落在局部最优上后两个模型的权重差了快 0.2。从那以后我每次跑 BMA 都强制走一遍多初始值扫描 对数似然曲线检查 分时段权重对比的流程再进结果解读。希望帮到你。本文还有配套的精品资源点击获取