简介这份资源围绕多模型集成中的BMABagging of Multiple Models方法展开面向具备一定MATLAB基础、从事气候预测、环境研究或统计建模的学习者与科研人员帮助解决单一模型预测不稳定、不确定性偏高的问题。压缩包共13个文件约5.23MB以8个m脚本为核心配合2个mat数据文件与2个txt说明文件覆盖模型拟合、权重计算、结果评估与可视化等环节。脚本中既有BMA核心算法实现也包含模型参数拟合、加权平均预测以及散点图、密度图等分析工具数据文件则提供预测与观测样本便于直接复现完整流程。已有839人学习下载说明该方向具有稳定的关注度。读者可借此理解BMA如何依据各模型在验证集上的表现分配权重并通过加权平均提升预测精度与稳定性同时掌握从数据准备、模型构建到结果评估的完整实现思路适合作为多模型集成方法的入门实践与二次开发参考。1. BMA 多模型集成把几个模型的预测揉成一个更稳的结果你手头有三个模型XGBoost 在结构化特征上表现不错LSTM 抓住了时间序列的局部模式一个线性回归作为兜底。单独跑AUC 忽高忽低换一批数据就翻车。你想把它们合起来用但简单平均权重全靠拍脑袋加权求和又找不到合理的权重依据。这时候 BMABayesian Model Averaging贝叶斯模型平均就派上用场了。BMA 的核心思路不复杂它不是选一个“最好”的模型而是根据每个模型在观测数据上的后验概率来加权把多个模型的预测分布合成一个混合分布。权重不是人为设定的而是从数据里推出来的。这意味着表现好的模型自动获得更高权重表现差的被压低而且当所有模型都烂的时候BMA 给出的预测不确定性会变大而不是硬给一个自信的错答案。这套方法适合谁做风控评分、销量预测、设备故障预警的工程师只要你有多个模型的预测输出又不想手动调集成权重BMA 就是一条值得走的路。下面从原理到代码把这条路铺清楚。2. BMA 的数学骨架与权重从哪来2.1 后验权重不是拍脑袋是算出来的BMA 的出发点就是贝叶斯定理。假设有 K 个候选模型 M₁ 到 M_K观测数据为 D那 BMA 对目标变量 y 的预测分布是P(y | D) Σₖ P(y | Mₖ, D) · P(Mₖ | D)其中 P(Mₖ | D) 就是模型 Mₖ 的后验权重也就是它在看到数据之后“应该占多大话语权”。这个后验权重又由两部分决定P(Mₖ | D) ∝ P(D | Mₖ) · P(Mₖ)P(D | Mₖ) 是边际似然也叫模型证据表示这个模型对数据的解释能力P(Mₖ) 是先验权重如果你对某个模型有额外信任可以调高没有就先设均匀分布。关键点在于边际似然 P(D | Mₖ) 的计算要对模型参数积分P(D | Mₖ) ∫ P(D | θₖ, Mₖ) · P(θₖ | Mₖ) dθₖ这个积分在高维参数空间里通常没有解析解所以实操中要么用 Laplace 近似、要么用 MCMC 采样、要么用 BIC 近似。我一般先用 BIC 近似快速跑通再决定要不要上 MCMC。2.2 用 BIC 近似边际似然最省事的起步方式BIC贝叶斯信息准则和边际似然之间有一个经典近似关系log P(D | Mₖ) ≈ -0.5 · BICₖ而 BICₖ -2 · log L̂ₖ pₖ · log(n)其中 L̂ₖ 是模型的最大似然值pₖ 是参数个数n 是样本量。把这个代入后验权重公式得到wₖ exp(-0.5 · BICₖ) / Σⱼ exp(-0.5 · BICⱼ)这就是实操中最常用的 BMA 权重计算方式。它不需要对参数做积分只需要每个模型的最大似然值和参数个数几乎任何模型都能给出这两个量。下面是一段可以直接跑的 Python 代码输入是每个模型的 log-likelihood 和参数数量输出 BMA 权重import numpy as np def bma_weights_bic(log_likelihoods, n_params, n_samples): 用 BIC 近似计算 BMA 后验权重。 参数: log_likelihoods: list[float]每个模型在验证集上的对数似然 n_params: list[int]每个模型的自由参数个数 n_samples: int样本量 返回: weights: np.ndarray归一化后的 BMA 权重 log_likelihoods np.array(log_likelihoods) n_params np.array(n_params) # 计算每个模型的 BIC bic -2 * log_likelihoods n_params * np.log(n_samples) # 用 log-sum-exp 技巧避免数值下溢 log_w -0.5 * bic log_w - np.max(log_w) # 减去最大值防止 exp 溢出 w np.exp(log_w) return w / w.sum() # 示例三个模型的 log-likelihood 和参数个数 ll [-120.5, -135.2, -128.7] params [12, 45, 8] n 500 weights bma_weights_bic(ll, params, n) print(BMA 权重:, weights)这段代码里有两个容易忽略的细节。第一log_w - np.max(log_w)这一步是数值稳定处理当 BIC 值很大时 exp 会下溢到 0减去最大值后至少有一个模型权重非零。第二log-likelihood 必须是在同一份验证集上算出来的不同数据集上的似然值没有可比性。参数说明log_likelihoods通常用验证集上的对数似然如果是分类任务用sklearn的log_loss乘以样本数取负即可n_params要算上所有可训练参数包括正则化项里的超参数不算在内n_samples是验证集样本量不是训练集。2.3 从权重到最终预测回归和分类的合成分开处理拿到权重之后合成预测的方式取决于任务类型。回归任务直接加权平均ŷ_BMA Σₖ wₖ · ŷₖ分类任务有两种做法。一种是加权平均概率再取 argmaxP̂(yc | x) Σₖ wₖ · P̂ₖ(yc | x)另一种是加权投票。我一般用概率平均因为它保留了不确定性信息后续做阈值调整更方便。def bma_predict_regression(predictions, weights): 回归任务的 BMA 合成预测。 参数: predictions: np.ndarray, shape(n_samples, n_models) weights: np.ndarray, shape(n_models,) 返回: y_pred: np.ndarray, shape(n_samples,) return predictions weights def bma_predict_classification(probas, weights): 分类任务的 BMA 合成概率。 参数: probas: np.ndarray, shape(n_samples, n_classes, n_models) weights: np.ndarray, shape(n_models,) 返回: proba_bma: np.ndarray, shape(n_samples, n_classes) # 对每个模型的概率按权重加权 proba_bma np.tensordot(probas, weights, axes([2], [0])) # 归一化防止浮点误差导致概率和不为 1 proba_bma / proba_bma.sum(axis1, keepdimsTrue) return proba_bmanp.tensordot那行是核心操作它把(n_samples, n_classes, n_models)和(n_models,)沿着模型维度做加权求和得到(n_samples, n_classes)。归一化那一步不能省因为浮点运算后概率和可能偏离 1影响后续 log_loss 计算。3. 把 BMA 集成跑通从基模型训练到权重验证3.1 基模型的选择与训练多样性比单模精度更重要BMA 的效果高度依赖基模型的多样性。如果三个模型都是 XGBoost 只是随机种子不同BMA 权重会集中在其中一个上集成收益几乎为零。我一般会选三类模型一个树模型XGBoost/LightGBM、一个线性模型Ridge/Logistic、一个神经网络或 KNN。它们在偏差-方差谱上处于不同位置BMA 才有得选。训练流程上每个基模型独立训练但必须用同一份训练集和验证集划分。我习惯用 5 折交叉验证每折训练一个基模型在验证折上记录预测值和 log-likelihood。这样每个模型有 5 组验证预测可以更稳健地估计权重。from sklearn.model_selection import KFold from sklearn.metrics import log_loss import numpy as np def train_base_models_cv(models, X, y, n_splits5): 用 K 折交叉验证训练基模型收集验证集预测和 log-likelihood。 参数: models: dict, {模型名: sklearn 模型实例} X: 特征矩阵 y: 标签 返回: cv_results: dict, 每个模型的验证预测和 log-likelihood kf KFold(n_splitsn_splits, shuffleTrue, random_state42) cv_results {name: {preds: [], loglik: []} for name in models} for train_idx, val_idx in kf.split(X): X_train, X_val X[train_idx], X[val_idx] y_train, y_val y[train_idx], y[val_idx] for name, model in models.items(): model_clone clone(model) model_clone.fit(X_train, y_train) if hasattr(model_clone, predict_proba): proba model_clone.predict_proba(X_val) ll -log_loss(y_val, proba, normalizeFalse) else: pred model_clone.predict(X_val) # 回归任务用高斯对数似然近似 resid y_val - pred sigma np.std(resid) 1e-8 ll -0.5 * np.sum((resid / sigma) ** 2 np.log(2 * np.pi * sigma ** 2)) cv_results[name][preds].append(proba if hasattr(model_clone, predict_proba) else pred) cv_results[name][loglik].append(ll) return cv_results这段代码里clone来自sklearn.base确保每折用的是未训练的新模型。log_loss的normalizeFalse返回的是总对数似然而非平均值这样不同折之间可以直接相加。回归任务用高斯似然近似sigma用残差标准差估计加1e-8防止除零。3.2 权重计算与集成预测的完整流程有了每折的 log-likelihood取平均后代入 BIC 公式算权重。然后对测试集做集成预测。from sklearn.base import clone def compute_bma_weights(cv_results, n_params_dict, n_samples): 从交叉验证结果计算 BMA 权重。 参数: cv_results: train_base_models_cv 的输出 n_params_dict: dict, {模型名: 参数个数} n_samples: 验证集样本量 返回: weights: dict, {模型名: BMA 权重} names list(cv_results.keys()) avg_ll [np.mean(cv_results[n][loglik]) for n in names] n_params [n_params_dict[n] for n in names] w bma_weights_bic(avg_ll, n_params, n_samples) return dict(zip(names, w)) def bma_ensemble_predict(models, weights, X_test, taskclassification): 训练最终基模型并做 BMA 集成预测。 参数: models: dict, {模型名: 模型实例} weights: dict, {模型名: BMA 权重} X_test: 测试集特征 task: classification 或 regression 返回: 集成预测结果 names list(models.keys()) w np.array([weights[n] for n in names]) if task classification: probas [] for name in names: model clone(models[name]) model.fit(X_train_full, y_train_full) probas.append(model.predict_proba(X_test)) probas np.stack(probas, axis2) # (n_samples, n_classes, n_models) return bma_predict_classification(probas, w) else: preds [] for name in names: model clone(models[name]) model.fit(X_train_full, y_train_full) preds.append(model.predict(X_test)) preds np.stack(preds, axis1) # (n_samples, n_models) return bma_predict_regression(preds, w)注意bma_ensemble_predict里用了X_train_full和y_train_full这是全量训练集。交叉验证阶段只是为了估权重最终模型要用全部数据重新训练。这个流程我踩过坑有一次直接用交叉验证的模型做预测结果每个折的模型只见过 80% 数据集成效果比单模还差。3.3 验证 BMA 是否真的比单模好别只看均值验证不能只看集成后的 AUC 或 RMSE 均值还要看稳定性和校准度。我一般做三件事第一对比 BMA 和最优单模在测试集上的指标如果提升不到 1%要怀疑基模型多样性不够。第二看 BMA 权重的分布如果某个模型权重超过 0.9说明其他模型没提供增量信息。第三画可靠性图reliability diagram检查 BMA 输出的概率是否校准。BMA 的一个理论优势是不确定性量化更合理如果校准反而变差说明权重估计有问题。from sklearn.calibration import calibration_curve import matplotlib.pyplot as plt def plot_calibration(y_true, proba_bma, proba_best_single, n_bins10): 对比 BMA 和最优单模的校准曲线。 fig, ax plt.subplots(figsize(6, 6)) for proba, label in [(proba_best_single, Best Single), (proba_bma, BMA)]: prob_true, prob_pred calibration_curve(y_true, proba[:, 1], n_binsn_bins) ax.plot(prob_pred, prob_true, markero, labellabel) ax.plot([0, 1], [0, 1], k--, labelPerfectly calibrated) ax.set_xlabel(Predicted probability) ax.set_ylabel(Fraction of positives) ax.legend() plt.tight_layout() plt.savefig(calibration_bma.png, dpi150)如果 BMA 的校准曲线比单模更贴近对角线说明集成确实改善了不确定性估计。如果更差检查权重计算时用的 log-likelihood 是否来自同一份数据以及基模型的概率输出是否已经过校准。4. BMA 集成避坑五个让我返工的血泪教训4.1 基模型 log-likelihood 不在同一数据集上算权重全乱现象BMA 权重算出来某个模型占 0.95但那个模型单独跑测试集指标并不好。原因不同模型的 log-likelihood 是在不同折或不同预处理后的数据上算的。比如模型 A 用了标准化后的特征模型 B 用了原始特征验证集样本量虽然一样但似然值的尺度不可比。解决所有基模型必须用完全相同的训练集、验证集划分和预处理流程。我现在的做法是把预处理放进 Pipeline确保每个模型看到的数值完全一致。log-likelihood 统一在验证折上计算不混用训练集似然。4.2 参数个数漏算BIC 惩罚项偏小现象复杂模型比如深层神经网络的 BMA 权重异常高但它在测试集上过拟合严重。原因n_params只算了权重数量漏了偏置项、批归一化参数、嵌入层参数。BIC 的惩罚项 p·log(n) 被低估复杂模型没被罚够。解决用sum(p.numel() for p in model.parameters())统计 PyTorch 模型参数量sklearn 模型用model.get_params()配合文档估算。宁可多算一点不要让复杂模型白占便宜。4.3 分类任务用 argmax 投票代替概率平均现象BMA 集成后的 AUC 比单模还低。原因每个基模型先 argmax 变成硬标签再投票。这样丢失了概率信息三个模型都预测类别 1 但概率分别是 0.51、0.52、0.53和概率 0.9、0.85、0.88 的投票结果一样但后者显然更可信。解决分类任务一律用概率平均保留每个模型的 softmax 输出。如果某个模型没有predict_proba用 Platt scaling 或 isotonic regression 先校准出概率。4.4 权重在训练集上算测试集上失效现象交叉验证时 BMA 比单模好 3%上线后反而差 2%。原因权重是用训练集上的 log-likelihood 算的训练集和线上数据分布有偏移训练集上表现好的模型在线上不一定好。解决权重计算必须用验证集或留出集不能用训练集。如果数据分布随时间变化用时间序列交叉验证最近的时间折权重更高。我一般会保留一个最近的验证集专门算权重不参与任何模型训练。4.5 基模型数量太多边际收益递减还引入噪声现象从 3 个模型加到 10 个模型BMA 效果没提升推理时间翻了 3 倍。原因新增的模型和已有模型高度相关BMA 权重被分散到相似模型上相当于降低了有效样本量。而且弱模型的 log-likelihood 低但权重不会为零仍然贡献噪声。解决基模型控制在 3 到 5 个选相关性低的。算一下模型间预测值的相关系数矩阵如果两个模型相关系数超过 0.95只保留一个。弱模型如果权重低于 0.05直接从集成里去掉重新归一化权重。5. 用 MCMC 采样替代 BIC 近似什么时候值得多花算力BIC 近似在样本量足够大、模型接近真实分布时表现不错但它有两个硬伤一是假设参数后验是高斯分布二是没法给出权重的置信区间。当你需要知道“这个权重到底有多确定”或者样本量小于 200 时MCMC 采样是更靠谱的选择。我一般用 PyMC 或 Stan 做 BMA 的 MCMC 实现。核心思路是把模型索引当成一个离散隐变量用 Gibbs 采样交替更新模型索引和模型参数。下面是一个用 PyMC 做 BMA 的简化示例假设有两个线性回归模型import pymc as pm import numpy as np # 模拟数据 np.random.seed(42) n 150 X1 np.random.randn(n, 2) X2 np.random.randn(n, 5) y 2.0 * X1[:, 0] - 1.5 * X1[:, 1] np.random.randn(n) * 0.5 with pm.Model() as bma_model: # 模型索引0 表示模型 11 表示模型 2 model_idx pm.Categorical(model_idx, p[0.5, 0.5]) # 模型 1 的参数 beta1 pm.Normal(beta1, mu0, sigma10, shape2) sigma1 pm.HalfNormal(sigma1, sigma1) mu1 pm.math.dot(X1, beta1) # 模型 2 的参数 beta2 pm.Normal(beta2, mu0, sigma10, shape5) sigma2 pm.HalfNormal(sigma2, sigma1) mu2 pm.math.dot(X2, beta2) # 根据模型索引选择均值 mu pm.math.switch(pm.math.eq(model_idx, 0), mu1, mu2) sigma pm.math.switch(pm.math.eq(model_idx, 0), sigma1, sigma2) # 似然 y_obs pm.Normal(y_obs, mumu, sigmasigma, observedy) # 采样 trace pm.sample(2000, tune1000, chains4, random_seed42) # 提取模型后验权重 model_idx_samples trace.posterior[model_idx].values.flatten() w_model1 np.mean(model_idx_samples 0) w_model2 np.mean(model_idx_samples 1) print(f模型 1 权重: {w_model1:.3f}, 模型 2 权重: {w_model2:.3f})这段代码的关键在pm.math.switch它根据model_idx的值动态选择不同模型的均值函数。采样完成后model_idx的后验分布直接给出了每个模型的权重而且还能看到权重的置信区间。MCMC 的代价是计算时间。上面这个例子在普通笔记本上跑 4 条链各 2000 次采样大约需要 2 到 3 分钟而 BIC 近似是毫秒级。我的经验是样本量大于 500 且模型参数少于 50 个时BIC 和 MCMC 的权重差异通常在 0.02 以内用 BIC 就够了样本量小或者需要权重的不确定性估计时才上 MCMC。还有一个实操技巧先用 BIC 算一遍权重把权重低于 0.01 的模型直接剔除再用 MCMC 对剩下的模型精细采样。这样既省时间又避免了弱模型干扰采样效率。最后说一个我自己的习惯每次做完 BMA 集成我都会把权重、基模型列表和验证集指标存成一个 JSON 文件和模型文件放在一起。过三个月回头看能清楚知道当时为什么选这几个模型、权重是多少、验证集表现如何。这个习惯帮我省了好几次“这个模型到底怎么来的”的返工。希望帮到你。本文还有配套的精品资源点击获取