1. 项目概述为什么我们需要CRPS在贝叶斯机器学习的实战中我们常常面临一个核心挑战如何评价一个模型的好坏这听起来简单但当你手上拿到的不是一个单一的预测值而是一个完整的概率分布时问题就变得复杂了。传统的点预测评估指标比如均方误差MSE或平均绝对误差MAE在这里就“失灵”了。它们只关心预测的“中心”离真实值有多远却完全忽略了模型对整个不确定性范围的刻画是否准确。一个预测分布很“瘦”过于自信但中心点偏得离谱的模型和一个预测分布很“胖”覆盖了各种可能性且中心点接近的模型用MSE去评价前者可能得分更高但这显然不是我们想要的。我们真正需要的是一个能同时评估预测分布的“位置”和“形状”的评分规则。这就是连续分级概率评分Continuous Ranked Probability Score, CRPS登场的背景。它不是一个象牙塔里的理论玩具而是数据科学家和机器学习工程师在评估贝叶斯回归、概率预测、气象预报乃至金融风险模型时手中一把不可或缺的“量尺”。简单说CRPS衡量的是你预测的整个概率分布与最终那个单一的现实观测值之间的“距离”距离越小说明你的概率预测越精准、越可靠。2. CRPS的核心原理与数学直觉要理解CRPS我们不能只停留在公式层面更需要建立起直观的“几何”感觉。它本质上是在比较两个函数你模型给出的预测累积分布函数CDF和观测值的“理想”CDF。2.1 从累积分布函数CDF说起任何一个概率分布无论是高斯分布、学生t分布还是复杂的混合模型都可以用它的累积分布函数 $F(y)$ 来完整描述。$F(y)$ 表示随机变量 $Y$ 取值小于等于 $y$ 的概率。它是一个从0单调递增到1的函数。 对于一次真实的观测 $x$一个具体的数值我们可以定义一个“理想”的CDF即海维赛德阶跃函数 $H(y-x)$ $$ H(y-x) \begin{cases} 0, \text{if } y x \ 1, \text{if } y \geq x \end{cases} $$ 这个函数在 $yx$ 处从0瞬间跳到1它完美地描述了“观测值 $x$ 以100%的概率小于等于任何大于等于 $x$ 的值而以0%的概率小于任何小于 $x$ 的值”。这是一个确定性事件的CDF。2.2 CRPS的定义一种“距离”CRPS就是衡量你的预测CDF $F$ 与这个理想阶跃CDF $H$ 之间差异的积分。其定义式为 $$ CRPS(F, x) \int_{-\infty}^{\infty} [F(y) - H(y-x)]^2 dy $$ 你可以把它想象成在整条实数轴上逐点计算两个函数值之差的平方然后求和积分。这个值永远是非负的只有当你的预测CDF $F$ 完美地与阶跃函数 $H$ 重合时CRPS才为0但这在概率预测中几乎不可能发生因为模型总会包含不确定性。为什么是平方积分这继承了平方误差的良好数学性质它是对称的、可微的并且对大的偏差给予更大的惩罚。在评估预测分布时这意味着如果模型预测的概率在真实值附近累积得太慢或太快都会导致 $F(y)$ 与 $H(y-x)$ 之间产生一个持续存在的差异区域从而被积分捕捉到形成较大的CRPS值。2.3 一个更易计算的等价形式直接计算上述积分通常不方便。一个更实用的等价形式是 $$ CRPS(F, x) E_F|Y - x| - \frac{1}{2} E_F|Y - Y| $$ 其中 $Y$ 和 $Y$ 是独立同分布都服从预测分布 $F$。这个公式提供了极其深刻的直觉第一项 $E_F|Y - x|$ 预测分布 $F$ 中随机样本 $Y$ 到观测值 $x$ 的平均绝对距离。这项衡量的是预测分布的“位置”准不准。如果分布的中心远离 $x$这项会很大。第二项 $-\frac{1}{2} E_F|Y - Y|$ 预测分布 $F$ 内部两个独立样本之间平均绝对距离的一半的负值。这项衡量的是预测分布本身的“离散度”或“宽度”。分布越分散越不确定两个样本间的平均距离就越大但前面有个负号所以减去一个更大的数会使最终的CRPS值变小。核心直觉CRPS鼓励模型在“位置准确”第一项小和“不确定性校准合理”第二项大之间取得最佳平衡。一个过于自信方差很小但位置有偏的模型第一项大第二项小因为样本间距离小CRPS大。一个过于保守方差很大但位置正确的模型第一项可能中等因为很多样本离 $x$ 远但第二项很大相减后CRPS可能反而较小。但如果方差大到离谱第一项也会被拉大。一个位置准确且不确定性适中的模型第一项小第二项适中CRPS最小。注意CRPS是“负向指标”即分数越低越好这与MSE、MAE等损失函数一致。3. 不同预测形式下的CRPS计算实战理论很美但落地到代码里才是关键。在实际的贝叶斯机器学习中我们得到的预测分布 $F$ 通常不是一个有解析表达式的简单分布而是以后验样本的形式存在。下面我们分场景讨论如何计算CRPS。3.1 场景一解析分布如高斯分布当预测分布是参数化分布时CRPS可能有闭合解。对于高斯分布 $N(\mu, \sigma^2)$其CRPS有一个相对简洁的表达式 $$ CRPS(N(\mu, \sigma^2), x) \sigma \left[ \frac{x-\mu}{\sigma} \left( 2\Phi\left(\frac{x-\mu}{\sigma}\right) - 1 \right) 2\phi\left(\frac{x-\mu}{\sigma}\right) - \frac{1}{\sqrt{\pi}} \right] $$ 其中 $\Phi$ 是标准正态CDF$\phi$ 是标准正态PDF。Python实现示例import numpy as np from scipy.stats import norm def crps_gaussian(mu, sigma, x): 计算高斯分布预测对于观测x的CRPS。 参数: mu: 预测均值 sigma: 预测标准差 x: 观测值 返回: CRPS值 # 标准化误差 z (x - mu) / sigma # 标准正态的PDF和CDF phi norm.pdf(z) Phi norm.cdf(z) # 应用公式 crps sigma * (z * (2 * Phi - 1) 2 * phi - 1 / np.sqrt(np.pi)) return crps # 示例 mu_pred, sigma_pred 5.0, 2.0 x_obs 4.2 print(fCRPS: {crps_gaussian(mu_pred, sigma_pred, x_obs):.4f})3.2 场景二经验分布后验样本这是贝叶斯机器学习中最常见的情况。我们通过MCMC如PyMC3、Stan或变分推断得到参数的后验样本进而得到预测分布的后验样本 ${y^{(1)}, y^{(2)}, ..., y^{(S)}}$其中 $S$ 是样本量。此时的预测CDF $F$ 是一个经验分布函数。对于经验分布最稳定高效的计算方法是利用其等价形式 $$ CRPS(F, x) \frac{1}{S} \sum_{s1}^{S} |y^{(s)} - x| - \frac{1}{2S^2} \sum_{s1}^{S} \sum_{t1}^{S} |y^{(s)} - y^{(t)}| $$ 这个公式直接翻译了 $E_F|Y-x|$ 和 $E_F|Y-Y|$ 的样本估计。但双重求和的计算复杂度是 $O(S^2)$当后验样本量 $S$ 很大时例如上万计算会非常慢。优化计算我们可以利用排序来将计算复杂度降至 $O(S \log S)$。对后验样本进行排序记为 $y_{(1)} \leq y_{(2)} \leq ... \leq y_{(S)}$则CRPS可以计算为 $$ CRPS \frac{1}{S} \sum_{s1}^{S} |y_{(s)} - x| - \frac{1}{S^2} \sum_{s1}^{S} (S - s 0.5) \cdot y_{(s)} $$ 注意这个公式需要一些推导更通用的做法是使用专门的库。Python实现与库推荐 对于经验分布强烈推荐使用properscoring库或xarray生态的xskillscore库它们经过优化且支持多维数组计算。import numpy as np import properscoring as ps # 生成后验预测样本 (S1000) np.random.seed(42) S 1000 # 假设预测分布是一个混合分布0.7*N(5,1) 0.3*N(10,2) component np.random.choice([0,1], sizeS, p[0.7, 0.3]) posterior_samples np.where(component 0, np.random.normal(5, 1, S), np.random.normal(10, 2, S)) x_obs 6.5 # 使用 properscoring 计算CRPS crps_val ps.crps_ensemble(x_obs, posterior_samples) print(fCRPS (using properscoring): {crps_val:.4f}) # 验证使用原始定义式慢仅用于小样本验证 def crps_empirical_naive(samples, x): S len(samples) term1 np.mean(np.abs(samples - x)) # 计算所有样本对之间的绝对差 diff_matrix np.abs(samples[:, None] - samples[None, :]) term2 0.5 * np.mean(diff_matrix) return term1 - term2 crps_naive crps_empirical_naive(posterior_samples, x_obs) print(fCRPS (naive double sum): {crps_naive:.4f})你会看到两个结果非常接近。properscoring内部使用了更高效的算法。3.3 场景三参数分布的集成多模型或混合后验有时我们会有多个候选模型或者后验本身是多个模式的混合。一种稳健的做法是计算每个模型/组件预测分布的CRPS然后取其平均集成CRPS。这等价于用集成分布作为最终的预测分布 $F$。 假设我们有 $K$ 个预测分布 $F_1, F_2, ..., F_K$集成分布为 $F_{avg} \frac{1}{K}\sum_{k1}^K F_k$那么 $$ CRPS(F_{avg}, x) \frac{1}{K} \sum_{k1}^{K} CRPS(F_k, x) - \frac{1}{2K^2} \sum_{k1}^{K}\sum_{l1}^{K} \int |F_k(y) - F_l(y)| dy $$ 第二项是各分布之间差异的惩罚项。在实践中如果各分布样本都可得最直接的方法是将所有 $K$ 个分布的样本合并成一个大样本集然后将其视为一个经验分布用3.2节的方法计算CRPS。4. CRPS在贝叶斯工作流中的应用与评估策略理解了如何计算单个预测-观测对的CRPS后我们需要将其应用到整个模型评估流程中。通常我们会在一个测试集 ${(x_i, y_i)}_{i1}^N$ 上评估模型。4.1 整体模型评估平均CRPS最常用的指标是测试集上的平均CRPS $$ \overline{CRPS} \frac{1}{N} \sum_{i1}^{N} CRPS(F_i, y_i) $$ 其中 $F_i$ 是模型在输入 $x_i$ 处给出的预测分布。$\overline{CRPS}$ 提供了一个模型整体概率预测校准性与锐度的综合度量。我们可以用它来比较不同贝叶斯模型例如不同先验、不同似然函数、不同推断算法的优劣。实操步骤模型训练与后验采样在训练集上训练你的贝叶斯模型如贝叶斯线性回归、高斯过程回归、贝叶斯神经网络并获得后验分布。生成后验预测分布对于测试集中的每一个样本 $x_i$利用后验分布生成其对应的预测分布 $F_i$。这通常意味着从后验中抽取参数样本然后为每个参数样本计算一个可能的 $y_i$ 值从而形成 $y_i$ 的后验预测样本集。逐点计算CRPS对每个测试点 $(x_i, y_i)$根据 $F_i$ 的形式解析或经验使用第3节的方法计算 $CRPS(F_i, y_i)$。计算平均对所有测试点的CRPS取平均。# 假设我们已经有了测试集 predictions_list 和 observations_list # predictions_list[i] 是对应于 observations_list[i] 的后验预测样本形状为 (S,) # observations_list[i] 是标量观测值 def evaluate_model_crps(predictions_list, observations_list): 计算模型在测试集上的平均CRPS。 total_crps 0.0 N len(observations_list) for i in range(N): crps_i ps.crps_ensemble(observations_list[i], predictions_list[i]) total_crps crps_i average_crps total_crps / N return average_crps # 示例假设有3个测试点 test_predictions [ np.random.normal(5, 1, 1000), # 对第一个点的预测样本 np.random.normal(8, 1.5, 1000), # 对第二个点 np.random.normal(6, 0.8, 1000) # 对第三个点 ] test_observations [5.2, 7.8, 6.5] avg_crps evaluate_model_crps(test_predictions, test_observations) print(fAverage CRPS on test set: {avg_crps:.4f})4.2 与其他评分规则的对比CRPS并非唯一的概率评分规则。理解它的相对位置有助于我们正确选用。评分规则核心思想优点缺点适用场景CRPS比较预测CDF与理想阶跃CDF的 $L^2$ 距离。1. 对位置和尺度都敏感。2. 与MAE同单位易于解释。3. 对经验分布有稳定估计算法。1. 计算比对数评分略复杂。2. 对分布尾部的极端事件相对不敏感。通用首选尤其适用于连续变量的概率预测评估如回归、气象预报。对数评分计算观测值在预测概率密度函数PDF下的对数似然。1. 具有严格适普性鼓励报告真实信念。2. 对分布的所有细节包括尾部都敏感。1. 对预测PDF的形态假设敏感若预测PDF在观测处为0或极小得分会趋于负无穷。2. 单位是“纳特/比特”解释性稍差。理论分析、模型比较当预测分布是解析形式且远离零概率区域时。间隔评分评估某个置信区间如90%是否包含了真实值并对区间宽度进行惩罚。1. 直观易于向非专业人士解释。2. 直接关注特定置信水平下的校准。1. 只利用了分布的局部信息特定分位数。2. 需要选择置信水平不同水平结果可能不同。需要报告特定置信区间可靠性的场景如金融风险在险价值VaR评估。选择建议对于一般性的贝叶斯回归模型评估CRPS通常是更稳健和实用的选择。它避免了对数评分中可能出现的数值问题如预测PDF为零并且其“平均绝对误差”的物理单位让业务方更容易理解。如果你特别关心模型对极端事件的预测能力如金融中的黑天鹅事件则应同时查看对数评分和CRPS因为对数评分对分布尾部的概率质量更敏感。如果评估目标是特定置信区间的覆盖概率则间隔评分是更直接的工具。4.3 可视化诊断PIT直方图与分位数图平均CRPS是一个标量汇总指标但为了深入诊断模型概率预测的校准情况我们需要可视化工具。概率积分变换PIT直方图 对于每个测试点 $(x_i, y_i)$计算 $u_i F_i(y_i)$即观测值在其预测CDF下的分位数。如果模型完美校准那么 ${u_i}$ 应该服从均匀分布 $U(0,1)$。绘制 $u_i$ 的直方图平坦的直方图表明概率预测校准良好。U形直方图表明预测分布过于集中过于自信观测值常落在分布的尾部。倒U形拱形直方图表明预测分布过于分散过于保守观测值常聚集在分布中心。倾斜的直方图表明预测分布存在系统性偏差。import matplotlib.pyplot as plt import numpy as np def plot_pit_histogram(predictions_list, observations_list): 绘制PIT直方图以诊断概率预测的校准性。 pit_values [] for i in range(len(observations_list)): # 计算观测值在预测样本中的分位数 # np.mean(preds obs) 给出了经验CDF在obs处的值 pit np.mean(predictions_list[i] observations_list[i]) pit_values.append(pit) pit_values np.array(pit_values) plt.figure(figsize(10, 5)) plt.hist(pit_values, bins20, edgecolorblack, alpha0.7, densityTrue) plt.axhline(y1.0, colorr, linestyle--, labelPerfect Calibration (Uniform)) plt.xlabel(Probability Integral Transform (PIT)) plt.ylabel(Density) plt.title(PIT Histogram for Probabilistic Forecast Calibration) plt.legend() plt.show() # 使用之前的测试数据示例 plot_pit_histogram(test_predictions, test_observations)分位数-分位数图QQ图 将观测值的实际分位数与在完美校准下期望的分位数进行对比。如果点大致落在对角线上则校准良好。这些可视化工具与CRPS等标量指标结合构成了评估贝叶斯模型概率预测性能的完整工具箱。5. 高级话题与常见问题排查在实际应用中你可能会遇到一些棘手的情况。这里分享一些经验和解决方案。5.1 处理高维输出与多变量CRPS当预测目标 $y$ 是多维向量时例如空间多个点的温度预测我们需要一个能评估联合预测分布好坏的指标。直接推广的CRPS即积分在 $R^d$ 上进行计算会变得异常复杂。 一种常见且实用的替代方案是计算能量评分Energy Score它是CRPS在多维情况下的自然推广。对于后验样本 ${y^{(s)}}_{s1}^S$ 和观测 $x$能量评分定义为 $$ ES(F, x) E_F ||Y - x|| - \frac{1}{2} E_F ||Y - Y|| $$ 其中 $||\cdot||$ 通常是欧几里得范数。这完全类比了CRPS的等价形式只是将绝对值替换为范数。计算时同样用样本均值来估计期望。def energy_score(samples, observation, beta1.0): 计算多变量预测的能量评分。 参数: samples: 形状为 (S, d) 的预测样本数组S是样本数d是维度。 observation: 形状为 (d,) 的观测向量。 beta: 范数的指数通常为1欧几里得距离或2。 返回: 能量评分越低越好。 S, d samples.shape # 第一项样本到观测的平均距离 term1 np.mean(np.linalg.norm(samples - observation, ordbeta, axis1)) # 第二项样本间平均距离的一半 # 注意这里计算所有样本对的距离复杂度O(S^2)。对于大S可能需要随机采样样本对。 distances [] for i in range(S): for j in range(i1, S): # 避免重复和自身比较 dist np.linalg.norm(samples[i] - samples[j], ordbeta) distances.append(dist) term2 0.5 * np.mean(distances) if distances else 0.0 # 也可以使用向量化计算但可能内存消耗大 # diff samples[:, None, :] - samples[None, :, :] # (S, S, d) # pairwise_dist np.linalg.norm(diff, ordbeta, axis2) # term2 0.5 * np.mean(pairwise_dist[np.triu_indices(S, k1)]) return term1 - term25.2 CRPS的分解校准项与锐度项与Brier评分类似CRPS也可以被分解这有助于我们更精细地理解模型误差的来源。CRPS可以近似分解为 $$ CRPS \approx \text{可靠性} \text{锐度} $$可靠性Calibration衡量预测概率与观测频率的一致性。理想情况下当你预测某事件发生的概率为p时它实际发生的频率也应该是p。可靠性差意味着模型“说谎”了。锐度Sharpness衡量预测分布的集中程度。在保证可靠性的前提下预测分布越集中锐度越高说明模型越自信提供的信息量越大这是好的。这种分解可以通过对测试集样本进行分箱根据预测的某个特征如均值或方差来实现。在实践中我们更多是通过PIT直方图来定性评估可靠性而锐度则可以通过预测分布的方差或分位数的平均宽度来量化。5.3 常见陷阱与排查技巧后验样本量不足计算经验分布的CRPS时如果后验样本量 $S$ 太小例如少于100CRPS的估计会不稳定且方差很大。建议确保用于评估的预测样本量足够大通常 $S \geq 1000$。在MCMC采样中要确保链已收敛且采样充分。预测分布有偏如果模型的平均CRPS很高首先检查PIT直方图。如果是倾斜的说明预测存在系统性偏差。排查检查模型特征工程是否充分先验选择是否不合理似然函数是否假设错误例如实际数据有异方差性却用了同方差的高斯似然。预测分布过度自信或过度保守PIT直方图呈U形或倒U形。排查U形过度自信模型方差可能被低估。考虑使用更灵活的似然如学生t分布代替高斯分布或引入随机效应来捕捉未被解释的变异。倒U形过度保守模型方差可能被高估。检查是否有不相关的特征引入了噪声或者先验分布是否对参数赋予了过大的不确定性。计算效率问题当测试集很大$N$很大且每个点的后验样本量 $S$ 也很大时逐点计算CRPS可能较慢。优化使用向量化操作和优化过的库如properscoring。对于只需要比较模型相对性能的场景可以考虑在测试集上对数据进行下采样或者使用更少的后验样本进行计算需确保估计依然稳定。如果使用自定义实现确保使用了基于排序的 $O(S \log S)$ 算法而不是 $O(S^2)$ 的双重循环。与业务指标对齐CRPS是一个统计上良好的评分规则但最终模型的价值要体现在业务上。建议在报告CRPS的同时也计算并报告一些业务相关的分位数指标。例如在需求预测中除了平均CRPS还可以报告在90%分位数预测下的平均绝对误差因为这直接关系到库存成本。CRPS作为一把衡量概率预测好坏的“尺子”其价值在于它将我们对于“好预测”的直觉——既要准又要诚实地表达不确定性——转化为了一个可严格计算和优化的数学量。将它熟练地融入你的贝叶斯建模工作流中能让你对模型的认知从“点估计的迷雾”迈入“概率分布的清晰之境”。 贝叶斯机器学习中CRPS：评估概率预测的核心指标与实战指南
1. 项目概述为什么我们需要CRPS在贝叶斯机器学习的实战中我们常常面临一个核心挑战如何评价一个模型的好坏这听起来简单但当你手上拿到的不是一个单一的预测值而是一个完整的概率分布时问题就变得复杂了。传统的点预测评估指标比如均方误差MSE或平均绝对误差MAE在这里就“失灵”了。它们只关心预测的“中心”离真实值有多远却完全忽略了模型对整个不确定性范围的刻画是否准确。一个预测分布很“瘦”过于自信但中心点偏得离谱的模型和一个预测分布很“胖”覆盖了各种可能性且中心点接近的模型用MSE去评价前者可能得分更高但这显然不是我们想要的。我们真正需要的是一个能同时评估预测分布的“位置”和“形状”的评分规则。这就是连续分级概率评分Continuous Ranked Probability Score, CRPS登场的背景。它不是一个象牙塔里的理论玩具而是数据科学家和机器学习工程师在评估贝叶斯回归、概率预测、气象预报乃至金融风险模型时手中一把不可或缺的“量尺”。简单说CRPS衡量的是你预测的整个概率分布与最终那个单一的现实观测值之间的“距离”距离越小说明你的概率预测越精准、越可靠。2. CRPS的核心原理与数学直觉要理解CRPS我们不能只停留在公式层面更需要建立起直观的“几何”感觉。它本质上是在比较两个函数你模型给出的预测累积分布函数CDF和观测值的“理想”CDF。2.1 从累积分布函数CDF说起任何一个概率分布无论是高斯分布、学生t分布还是复杂的混合模型都可以用它的累积分布函数 $F(y)$ 来完整描述。$F(y)$ 表示随机变量 $Y$ 取值小于等于 $y$ 的概率。它是一个从0单调递增到1的函数。 对于一次真实的观测 $x$一个具体的数值我们可以定义一个“理想”的CDF即海维赛德阶跃函数 $H(y-x)$ $$ H(y-x) \begin{cases} 0, \text{if } y x \ 1, \text{if } y \geq x \end{cases} $$ 这个函数在 $yx$ 处从0瞬间跳到1它完美地描述了“观测值 $x$ 以100%的概率小于等于任何大于等于 $x$ 的值而以0%的概率小于任何小于 $x$ 的值”。这是一个确定性事件的CDF。2.2 CRPS的定义一种“距离”CRPS就是衡量你的预测CDF $F$ 与这个理想阶跃CDF $H$ 之间差异的积分。其定义式为 $$ CRPS(F, x) \int_{-\infty}^{\infty} [F(y) - H(y-x)]^2 dy $$ 你可以把它想象成在整条实数轴上逐点计算两个函数值之差的平方然后求和积分。这个值永远是非负的只有当你的预测CDF $F$ 完美地与阶跃函数 $H$ 重合时CRPS才为0但这在概率预测中几乎不可能发生因为模型总会包含不确定性。为什么是平方积分这继承了平方误差的良好数学性质它是对称的、可微的并且对大的偏差给予更大的惩罚。在评估预测分布时这意味着如果模型预测的概率在真实值附近累积得太慢或太快都会导致 $F(y)$ 与 $H(y-x)$ 之间产生一个持续存在的差异区域从而被积分捕捉到形成较大的CRPS值。2.3 一个更易计算的等价形式直接计算上述积分通常不方便。一个更实用的等价形式是 $$ CRPS(F, x) E_F|Y - x| - \frac{1}{2} E_F|Y - Y| $$ 其中 $Y$ 和 $Y$ 是独立同分布都服从预测分布 $F$。这个公式提供了极其深刻的直觉第一项 $E_F|Y - x|$ 预测分布 $F$ 中随机样本 $Y$ 到观测值 $x$ 的平均绝对距离。这项衡量的是预测分布的“位置”准不准。如果分布的中心远离 $x$这项会很大。第二项 $-\frac{1}{2} E_F|Y - Y|$ 预测分布 $F$ 内部两个独立样本之间平均绝对距离的一半的负值。这项衡量的是预测分布本身的“离散度”或“宽度”。分布越分散越不确定两个样本间的平均距离就越大但前面有个负号所以减去一个更大的数会使最终的CRPS值变小。核心直觉CRPS鼓励模型在“位置准确”第一项小和“不确定性校准合理”第二项大之间取得最佳平衡。一个过于自信方差很小但位置有偏的模型第一项大第二项小因为样本间距离小CRPS大。一个过于保守方差很大但位置正确的模型第一项可能中等因为很多样本离 $x$ 远但第二项很大相减后CRPS可能反而较小。但如果方差大到离谱第一项也会被拉大。一个位置准确且不确定性适中的模型第一项小第二项适中CRPS最小。注意CRPS是“负向指标”即分数越低越好这与MSE、MAE等损失函数一致。3. 不同预测形式下的CRPS计算实战理论很美但落地到代码里才是关键。在实际的贝叶斯机器学习中我们得到的预测分布 $F$ 通常不是一个有解析表达式的简单分布而是以后验样本的形式存在。下面我们分场景讨论如何计算CRPS。3.1 场景一解析分布如高斯分布当预测分布是参数化分布时CRPS可能有闭合解。对于高斯分布 $N(\mu, \sigma^2)$其CRPS有一个相对简洁的表达式 $$ CRPS(N(\mu, \sigma^2), x) \sigma \left[ \frac{x-\mu}{\sigma} \left( 2\Phi\left(\frac{x-\mu}{\sigma}\right) - 1 \right) 2\phi\left(\frac{x-\mu}{\sigma}\right) - \frac{1}{\sqrt{\pi}} \right] $$ 其中 $\Phi$ 是标准正态CDF$\phi$ 是标准正态PDF。Python实现示例import numpy as np from scipy.stats import norm def crps_gaussian(mu, sigma, x): 计算高斯分布预测对于观测x的CRPS。 参数: mu: 预测均值 sigma: 预测标准差 x: 观测值 返回: CRPS值 # 标准化误差 z (x - mu) / sigma # 标准正态的PDF和CDF phi norm.pdf(z) Phi norm.cdf(z) # 应用公式 crps sigma * (z * (2 * Phi - 1) 2 * phi - 1 / np.sqrt(np.pi)) return crps # 示例 mu_pred, sigma_pred 5.0, 2.0 x_obs 4.2 print(fCRPS: {crps_gaussian(mu_pred, sigma_pred, x_obs):.4f})3.2 场景二经验分布后验样本这是贝叶斯机器学习中最常见的情况。我们通过MCMC如PyMC3、Stan或变分推断得到参数的后验样本进而得到预测分布的后验样本 ${y^{(1)}, y^{(2)}, ..., y^{(S)}}$其中 $S$ 是样本量。此时的预测CDF $F$ 是一个经验分布函数。对于经验分布最稳定高效的计算方法是利用其等价形式 $$ CRPS(F, x) \frac{1}{S} \sum_{s1}^{S} |y^{(s)} - x| - \frac{1}{2S^2} \sum_{s1}^{S} \sum_{t1}^{S} |y^{(s)} - y^{(t)}| $$ 这个公式直接翻译了 $E_F|Y-x|$ 和 $E_F|Y-Y|$ 的样本估计。但双重求和的计算复杂度是 $O(S^2)$当后验样本量 $S$ 很大时例如上万计算会非常慢。优化计算我们可以利用排序来将计算复杂度降至 $O(S \log S)$。对后验样本进行排序记为 $y_{(1)} \leq y_{(2)} \leq ... \leq y_{(S)}$则CRPS可以计算为 $$ CRPS \frac{1}{S} \sum_{s1}^{S} |y_{(s)} - x| - \frac{1}{S^2} \sum_{s1}^{S} (S - s 0.5) \cdot y_{(s)} $$ 注意这个公式需要一些推导更通用的做法是使用专门的库。Python实现与库推荐 对于经验分布强烈推荐使用properscoring库或xarray生态的xskillscore库它们经过优化且支持多维数组计算。import numpy as np import properscoring as ps # 生成后验预测样本 (S1000) np.random.seed(42) S 1000 # 假设预测分布是一个混合分布0.7*N(5,1) 0.3*N(10,2) component np.random.choice([0,1], sizeS, p[0.7, 0.3]) posterior_samples np.where(component 0, np.random.normal(5, 1, S), np.random.normal(10, 2, S)) x_obs 6.5 # 使用 properscoring 计算CRPS crps_val ps.crps_ensemble(x_obs, posterior_samples) print(fCRPS (using properscoring): {crps_val:.4f}) # 验证使用原始定义式慢仅用于小样本验证 def crps_empirical_naive(samples, x): S len(samples) term1 np.mean(np.abs(samples - x)) # 计算所有样本对之间的绝对差 diff_matrix np.abs(samples[:, None] - samples[None, :]) term2 0.5 * np.mean(diff_matrix) return term1 - term2 crps_naive crps_empirical_naive(posterior_samples, x_obs) print(fCRPS (naive double sum): {crps_naive:.4f})你会看到两个结果非常接近。properscoring内部使用了更高效的算法。3.3 场景三参数分布的集成多模型或混合后验有时我们会有多个候选模型或者后验本身是多个模式的混合。一种稳健的做法是计算每个模型/组件预测分布的CRPS然后取其平均集成CRPS。这等价于用集成分布作为最终的预测分布 $F$。 假设我们有 $K$ 个预测分布 $F_1, F_2, ..., F_K$集成分布为 $F_{avg} \frac{1}{K}\sum_{k1}^K F_k$那么 $$ CRPS(F_{avg}, x) \frac{1}{K} \sum_{k1}^{K} CRPS(F_k, x) - \frac{1}{2K^2} \sum_{k1}^{K}\sum_{l1}^{K} \int |F_k(y) - F_l(y)| dy $$ 第二项是各分布之间差异的惩罚项。在实践中如果各分布样本都可得最直接的方法是将所有 $K$ 个分布的样本合并成一个大样本集然后将其视为一个经验分布用3.2节的方法计算CRPS。4. CRPS在贝叶斯工作流中的应用与评估策略理解了如何计算单个预测-观测对的CRPS后我们需要将其应用到整个模型评估流程中。通常我们会在一个测试集 ${(x_i, y_i)}_{i1}^N$ 上评估模型。4.1 整体模型评估平均CRPS最常用的指标是测试集上的平均CRPS $$ \overline{CRPS} \frac{1}{N} \sum_{i1}^{N} CRPS(F_i, y_i) $$ 其中 $F_i$ 是模型在输入 $x_i$ 处给出的预测分布。$\overline{CRPS}$ 提供了一个模型整体概率预测校准性与锐度的综合度量。我们可以用它来比较不同贝叶斯模型例如不同先验、不同似然函数、不同推断算法的优劣。实操步骤模型训练与后验采样在训练集上训练你的贝叶斯模型如贝叶斯线性回归、高斯过程回归、贝叶斯神经网络并获得后验分布。生成后验预测分布对于测试集中的每一个样本 $x_i$利用后验分布生成其对应的预测分布 $F_i$。这通常意味着从后验中抽取参数样本然后为每个参数样本计算一个可能的 $y_i$ 值从而形成 $y_i$ 的后验预测样本集。逐点计算CRPS对每个测试点 $(x_i, y_i)$根据 $F_i$ 的形式解析或经验使用第3节的方法计算 $CRPS(F_i, y_i)$。计算平均对所有测试点的CRPS取平均。# 假设我们已经有了测试集 predictions_list 和 observations_list # predictions_list[i] 是对应于 observations_list[i] 的后验预测样本形状为 (S,) # observations_list[i] 是标量观测值 def evaluate_model_crps(predictions_list, observations_list): 计算模型在测试集上的平均CRPS。 total_crps 0.0 N len(observations_list) for i in range(N): crps_i ps.crps_ensemble(observations_list[i], predictions_list[i]) total_crps crps_i average_crps total_crps / N return average_crps # 示例假设有3个测试点 test_predictions [ np.random.normal(5, 1, 1000), # 对第一个点的预测样本 np.random.normal(8, 1.5, 1000), # 对第二个点 np.random.normal(6, 0.8, 1000) # 对第三个点 ] test_observations [5.2, 7.8, 6.5] avg_crps evaluate_model_crps(test_predictions, test_observations) print(fAverage CRPS on test set: {avg_crps:.4f})4.2 与其他评分规则的对比CRPS并非唯一的概率评分规则。理解它的相对位置有助于我们正确选用。评分规则核心思想优点缺点适用场景CRPS比较预测CDF与理想阶跃CDF的 $L^2$ 距离。1. 对位置和尺度都敏感。2. 与MAE同单位易于解释。3. 对经验分布有稳定估计算法。1. 计算比对数评分略复杂。2. 对分布尾部的极端事件相对不敏感。通用首选尤其适用于连续变量的概率预测评估如回归、气象预报。对数评分计算观测值在预测概率密度函数PDF下的对数似然。1. 具有严格适普性鼓励报告真实信念。2. 对分布的所有细节包括尾部都敏感。1. 对预测PDF的形态假设敏感若预测PDF在观测处为0或极小得分会趋于负无穷。2. 单位是“纳特/比特”解释性稍差。理论分析、模型比较当预测分布是解析形式且远离零概率区域时。间隔评分评估某个置信区间如90%是否包含了真实值并对区间宽度进行惩罚。1. 直观易于向非专业人士解释。2. 直接关注特定置信水平下的校准。1. 只利用了分布的局部信息特定分位数。2. 需要选择置信水平不同水平结果可能不同。需要报告特定置信区间可靠性的场景如金融风险在险价值VaR评估。选择建议对于一般性的贝叶斯回归模型评估CRPS通常是更稳健和实用的选择。它避免了对数评分中可能出现的数值问题如预测PDF为零并且其“平均绝对误差”的物理单位让业务方更容易理解。如果你特别关心模型对极端事件的预测能力如金融中的黑天鹅事件则应同时查看对数评分和CRPS因为对数评分对分布尾部的概率质量更敏感。如果评估目标是特定置信区间的覆盖概率则间隔评分是更直接的工具。4.3 可视化诊断PIT直方图与分位数图平均CRPS是一个标量汇总指标但为了深入诊断模型概率预测的校准情况我们需要可视化工具。概率积分变换PIT直方图 对于每个测试点 $(x_i, y_i)$计算 $u_i F_i(y_i)$即观测值在其预测CDF下的分位数。如果模型完美校准那么 ${u_i}$ 应该服从均匀分布 $U(0,1)$。绘制 $u_i$ 的直方图平坦的直方图表明概率预测校准良好。U形直方图表明预测分布过于集中过于自信观测值常落在分布的尾部。倒U形拱形直方图表明预测分布过于分散过于保守观测值常聚集在分布中心。倾斜的直方图表明预测分布存在系统性偏差。import matplotlib.pyplot as plt import numpy as np def plot_pit_histogram(predictions_list, observations_list): 绘制PIT直方图以诊断概率预测的校准性。 pit_values [] for i in range(len(observations_list)): # 计算观测值在预测样本中的分位数 # np.mean(preds obs) 给出了经验CDF在obs处的值 pit np.mean(predictions_list[i] observations_list[i]) pit_values.append(pit) pit_values np.array(pit_values) plt.figure(figsize(10, 5)) plt.hist(pit_values, bins20, edgecolorblack, alpha0.7, densityTrue) plt.axhline(y1.0, colorr, linestyle--, labelPerfect Calibration (Uniform)) plt.xlabel(Probability Integral Transform (PIT)) plt.ylabel(Density) plt.title(PIT Histogram for Probabilistic Forecast Calibration) plt.legend() plt.show() # 使用之前的测试数据示例 plot_pit_histogram(test_predictions, test_observations)分位数-分位数图QQ图 将观测值的实际分位数与在完美校准下期望的分位数进行对比。如果点大致落在对角线上则校准良好。这些可视化工具与CRPS等标量指标结合构成了评估贝叶斯模型概率预测性能的完整工具箱。5. 高级话题与常见问题排查在实际应用中你可能会遇到一些棘手的情况。这里分享一些经验和解决方案。5.1 处理高维输出与多变量CRPS当预测目标 $y$ 是多维向量时例如空间多个点的温度预测我们需要一个能评估联合预测分布好坏的指标。直接推广的CRPS即积分在 $R^d$ 上进行计算会变得异常复杂。 一种常见且实用的替代方案是计算能量评分Energy Score它是CRPS在多维情况下的自然推广。对于后验样本 ${y^{(s)}}_{s1}^S$ 和观测 $x$能量评分定义为 $$ ES(F, x) E_F ||Y - x|| - \frac{1}{2} E_F ||Y - Y|| $$ 其中 $||\cdot||$ 通常是欧几里得范数。这完全类比了CRPS的等价形式只是将绝对值替换为范数。计算时同样用样本均值来估计期望。def energy_score(samples, observation, beta1.0): 计算多变量预测的能量评分。 参数: samples: 形状为 (S, d) 的预测样本数组S是样本数d是维度。 observation: 形状为 (d,) 的观测向量。 beta: 范数的指数通常为1欧几里得距离或2。 返回: 能量评分越低越好。 S, d samples.shape # 第一项样本到观测的平均距离 term1 np.mean(np.linalg.norm(samples - observation, ordbeta, axis1)) # 第二项样本间平均距离的一半 # 注意这里计算所有样本对的距离复杂度O(S^2)。对于大S可能需要随机采样样本对。 distances [] for i in range(S): for j in range(i1, S): # 避免重复和自身比较 dist np.linalg.norm(samples[i] - samples[j], ordbeta) distances.append(dist) term2 0.5 * np.mean(distances) if distances else 0.0 # 也可以使用向量化计算但可能内存消耗大 # diff samples[:, None, :] - samples[None, :, :] # (S, S, d) # pairwise_dist np.linalg.norm(diff, ordbeta, axis2) # term2 0.5 * np.mean(pairwise_dist[np.triu_indices(S, k1)]) return term1 - term25.2 CRPS的分解校准项与锐度项与Brier评分类似CRPS也可以被分解这有助于我们更精细地理解模型误差的来源。CRPS可以近似分解为 $$ CRPS \approx \text{可靠性} \text{锐度} $$可靠性Calibration衡量预测概率与观测频率的一致性。理想情况下当你预测某事件发生的概率为p时它实际发生的频率也应该是p。可靠性差意味着模型“说谎”了。锐度Sharpness衡量预测分布的集中程度。在保证可靠性的前提下预测分布越集中锐度越高说明模型越自信提供的信息量越大这是好的。这种分解可以通过对测试集样本进行分箱根据预测的某个特征如均值或方差来实现。在实践中我们更多是通过PIT直方图来定性评估可靠性而锐度则可以通过预测分布的方差或分位数的平均宽度来量化。5.3 常见陷阱与排查技巧后验样本量不足计算经验分布的CRPS时如果后验样本量 $S$ 太小例如少于100CRPS的估计会不稳定且方差很大。建议确保用于评估的预测样本量足够大通常 $S \geq 1000$。在MCMC采样中要确保链已收敛且采样充分。预测分布有偏如果模型的平均CRPS很高首先检查PIT直方图。如果是倾斜的说明预测存在系统性偏差。排查检查模型特征工程是否充分先验选择是否不合理似然函数是否假设错误例如实际数据有异方差性却用了同方差的高斯似然。预测分布过度自信或过度保守PIT直方图呈U形或倒U形。排查U形过度自信模型方差可能被低估。考虑使用更灵活的似然如学生t分布代替高斯分布或引入随机效应来捕捉未被解释的变异。倒U形过度保守模型方差可能被高估。检查是否有不相关的特征引入了噪声或者先验分布是否对参数赋予了过大的不确定性。计算效率问题当测试集很大$N$很大且每个点的后验样本量 $S$ 也很大时逐点计算CRPS可能较慢。优化使用向量化操作和优化过的库如properscoring。对于只需要比较模型相对性能的场景可以考虑在测试集上对数据进行下采样或者使用更少的后验样本进行计算需确保估计依然稳定。如果使用自定义实现确保使用了基于排序的 $O(S \log S)$ 算法而不是 $O(S^2)$ 的双重循环。与业务指标对齐CRPS是一个统计上良好的评分规则但最终模型的价值要体现在业务上。建议在报告CRPS的同时也计算并报告一些业务相关的分位数指标。例如在需求预测中除了平均CRPS还可以报告在90%分位数预测下的平均绝对误差因为这直接关系到库存成本。CRPS作为一把衡量概率预测好坏的“尺子”其价值在于它将我们对于“好预测”的直觉——既要准又要诚实地表达不确定性——转化为了一个可严格计算和优化的数学量。将它熟练地融入你的贝叶斯建模工作流中能让你对模型的认知从“点估计的迷雾”迈入“概率分布的清晰之境”。