github：
https://github.com/stanfordmlgroup/ngboost

文档：
https://stanfordmlgroup.github.io/ngboost/3-interpretation.html

本篇文章A new gradient boosting method from Standford Research: NGBoost (Python code included)适合希望了解概率回归的新手，亮点在于NGBoost能够输出完整的概率分布，而不仅仅是点估计。这种方法利用自然梯度提升了模型的稳定性和训练效率，适合需要不确定性量化的场景。

---

---

在医疗保健、金融和天气预报等许多关键应用中，机器学习模型不仅需要预测单个数值，还需要量化不确定性。自然梯度提升（NGBoost） 是一种概率预测算法，它扩展了传统的梯度提升技术，以生成完整的预测分布，而不仅仅是点估计。通过利用_自然梯度_，NGBoost 能够更稳定、更高效地训练模型，这些模型不仅输出_均值_，还输出所选概率分布的_所有参数_，从而提供稳健的不确定性估计和更好的可解释性。

完整文章： [1910.03225] NGBoost: 自然梯度提升用于概率预测
引用： @misc{duan2020ngboostnaturalgradientboosting,
title={NGBoost: Natural Gradient Boosting for Probabilistic Prediction},
author={Tony Duan and Anand Avati and Daisy Yi Ding and Khanh K. Thai and Sanjay Basu and Andrew Y. Ng and Alejandro Schuler},year={2020}, eprint={1910.03225}, archivePrefix={arXiv},primaryClass={cs.LG},url={https://arxiv.org/abs/1910.03225}, }

1. 为什么选择概率回归？

大多数回归模型产生的是点估计——给定输入特征的目标变量的平均值或期望值。然而，许多现实世界的问题需要更丰富、更全面的输出：对潜在结果的_完整概率分布_。这种分布预测支持：

• 不确定性估计（例如，置信区间），
• 风险评估（例如，预测极端结果），
• 不确定性下的决策制定（例如，预期损失、概率阈值）。

例如，在医学预后中，知道患者在未来六个月内有 70% 的生存几率，远比简单地预测他们能活 180 天更有用。

2. 现有解决方案及其局限性

有几种方法尝试进行概率回归：

• 贝叶斯模型（例如，BART、贝叶斯神经网络）自然地产生不确定性估计，但需要计算成本高昂的采样方法和专业知识才能有效实现。
• GAMLSS 允许异方差建模，但结构僵化。
• 点估计回归器 附加事后方差估计，依赖于同方差性等强假设，这在真实数据中很少成立。

这些方法通常缺乏模块化、可扩展性或易用性。NGBoost 直接解决了这些挑战。

3. 什么是 NGBoost？

NGBoost（自然梯度提升）是一种模块化、灵活的概率预测算法。它结合了梯度提升机（GBM） 的优势和自然梯度的严谨性，使模型能够输出完整的概率分布——而不仅仅是单一预测。

NGBoost 的核心在于将条件分布的参数（例如，正态分布的均值和方差）视为多参数提升的目标。它通过使用自然梯度——一种几何感知的优化方法——来改进传统的梯度提升，以增强收敛性和稳定性。

NGBoost 的关键组成部分：

• 基学习器：可以是任何可微分模型（通常是决策树）。
• 分布族：正态分布、拉普拉斯分布以及其他具有连续参数的分布。
• 评分规则：最大似然估计（MLE）、连续排名概率分数（CRPS）等。

4. 技术基础

4.1 适当评分规则

在概率回归中，目标是最小化评分规则而不是损失函数。适当评分规则为更好的概率预测分配较低的分数。

• 对数评分（MLE）：
  $L(y)=-\log P(y)$
• CRPS（连续排名概率分数）：
  一种评估完整累积分布函数（CDF）的稳健替代方法。

4.2 自然梯度

标准梯度下降对参数化敏感，在优化分布参数时可能导致低效的更新。NGBoost 使用自然梯度，它根据参数空间的基础几何形状进行调整。这使得优化：

• 对重参数化不变
• 更符合分布之间的实际“距离”

从数学上讲，自然梯度考虑了 Fisher 信息矩阵，该矩阵捕获了参数空间的曲率。

5. NGBoost 贡献总结

1. 多参数提升：使用独立的梯度提升更新，同时建模分布的多个参数（例如，均值和方差）。
2. 自然梯度：改进收敛动态，并纠正由朴素参数更新引起的失真。
3. 模块化：用户可以指定基学习器、分布类型和评分规则以适应其应用。
4. 经验性能：NGBoost 在各种领域中，在点估计和不确定性估计方面都表现出竞争力。

6. NGBoost 的实际应用

凭借其灵活性和即插即用设计，NGBoost 具有广泛的适用性。用例包括：

• 医学预测：预测事件发生时间（例如，生存分析）。
• 天气预报：输出未来温度的概率区间。
• 金融建模：估计风险和收益分布。

NGBoost 的一个主要优势是它能够很好地扩展到大型数据集，并且易于与 scikit-learn 风格的管道集成。它也不需要大量的超参数调优或概率建模方面的领域特定专业知识。

NGBoost 是任何需要不确定性感知预测任务的强大工具。它将梯度提升的多功能性与通过自然梯度进行的概率建模的原则性方法相结合。随着对可解释、可靠机器学习的需求增长——特别是在高风险领域——NGBoost 提供了一个引人注目的解决方案。

7 代码示例

7.1 通过 pip 安装：

pip install --upgrade ngboost

7.2 简单的 XGBoost 风格用法：

import numpy as np
import pandas as pd
from ngboost import NGBRegressor
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error


data_url = "http://lib.stat.cmu.edu/datasets/boston"
raw_df = pd.read_csv(data_url, sep=r"\s+", skiprows=22, header=None)

X = np.hstack([raw_df.values[::2, :], raw_df.values[1::2, :2]])
Y = raw_df.values[1::2, 2]

X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=42)

ngb = NGBRegressor()
ngb.fit(X_train, Y_train)

Y_preds = ngb.predict(X_test)
Y_dists = ngb.pred_dist(X_test)

test_MSE = mean_squared_error(Y_test, Y_preds)
print("Test MSE:", test_MSE)
test_NLL = -Y_dists.logpdf(Y_test).mean()
print("Test NLL:", test_NLL)

8 模型解释

来自文档：
https://stanfordmlgroup.github.io/ngboost/3-interpretation.html

获得特征重要性

ngb = NGBRegressor(verbose=False).fit(X_reg_train, Y_reg_train)

## Feature importance for loc trees
feature_importance_loc = ngb.feature_importances_[0]

## Feature importance for scale trees
feature_importance_scale = ngb.feature_importances_[1]

import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns

df_loc = pd.DataFrame({'feature':load_boston()['feature_names'],
                       'importance':feature_importance_loc})\
    .sort_values('importance',ascending=False)
df_scale = pd.DataFrame({'feature':load_boston()['feature_names'],
                       'importance':feature_importance_scale})\
    .sort_values('importance',ascending=False)

fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(13,6))
fig.suptitle("Feature importance plot for distribution parameters", fontsize=17)
sns.barplot(x='importance',y='feature',ax=ax1,data=df_loc, color="skyblue").set_title('loc param')
sns.barplot(x='importance',y='feature',ax=ax2,data=df_scale, color="skyblue").set_title('scale param')

使用SHAP 进行模型解释

import shap
shap.initjs()

## SHAP plot for loc trees
explainer = shap.TreeExplainer(ngb, model_output=0) # use model_output = 1 for scale trees
shap_values = explainer.shap_values(X_reg_train)
shap.summary_plot(shap_values, X_reg_train, feature_names=load_boston()['feature_names'])

9 调优

9.1 分阶段预测

所有的 NGBoost 对象都支持阶段预测。

ngb_cls = NGBClassifier(Dist=k_categorical(3), Score=LogScore, n_estimators=500, verbose=False).fit(X_cls_train, Y_cls_train)

例如，要在拟合了 415 个基学习器后获取前 5 个示例的预测，可以使用：

preds = ngb_cls.staged_predict(X_cls_test)
preds[415][0:5]

pred_dists = ngb_cls.staged_pred_dist(X_cls_test)
pred_dists[415][0:5].params

{'p0': array([0.99074995, 0.91368635, 0.00517919, 0.00517919, 0.00517919]),
 'p1': array([0.00860966, 0.03267806, 0.99450359, 0.99450359, 0.99450359]),
 'p2': array([0.00064039, 0.05363559, 0.00031722, 0.00031722, 0.00031722])}

这在结合验证集上的错误追踪时非常有用，可以通过传递 X_val 和 Y_val 参数来实现，然后检查 .best_val_loss_itr 实例属性。

ngb = NGBRegressor()
ngb.fit(X_reg_train, Y_reg_train, X_val=X_reg_test, Y_val=Y_reg_test)  # 在自己的工作中使用验证集而不是测试集
print(ngb.best_val_loss_itr)
best_preds = ngb.predict(X_reg_test, max_iter=ngb.best_val_loss_itr)

[iter 0] loss=3.6556 val_loss=3.5575 scale=0.5000 norm=3.4142
[iter 100] loss=3.1118 val_loss=3.1284 scale=1.0000 norm=3.9174
[iter 200] loss=2.4839 val_loss=2.6398 scale=2.0000 norm=4.0907
[iter 300] loss=2.0183 val_loss=2.7162 scale=1.0000 norm=1.5637
[iter 400] loss=1.8111 val_loss=3.1315 scale=1.0000 norm=1.3983

9.2 Early Stopping

NGBoost 也支持早期停止。如果在 fit() 中传递一个整数 early_stopping_rounds 和验证集 (X_val, Y_val)，算法将在验证损失连续 early_stopping_rounds 次迭代增加后停止运行。

_ = NGBRegressor().fit(X_reg_train, Y_reg_train, X_val=X_reg_test, Y_val=Y_reg_test, early_stopping_rounds=2)

[iter 0] loss=3.6556 val_loss=3.5575 scale=0.5000 norm=3.4142
[iter 100] loss=3.1118 val_loss=3.1292 scale=1.0000 norm=3.9174
[iter 200] loss=2.4839 val_loss=2.6422 scale=2.0000 norm=4.0907
== Early stopping achieved.
== Best iteration / VAL234 (val_loss=2.5693)

可以通过 val_sample_weight 参数传递验证集样本权重以供 fit 使用。

9.3 使用 sklearn 模型选择

sksklearn 方法与 NGBoost 兼容。

from sklearn.model_selection import GridSearchCV
from sklearn.tree import DecisionTreeRegressor

b1 = DecisionTreeRegressor(criterion='friedman_mse', max_depth=2)
b2 = DecisionTreeRegressor(criterion='friedman_mse', max_depth=4)

param_grid = {
    'minibatch_frac': [1.0, 0.5],
    'Base': [b1, b2]
}

ngb = NGBRegressor(Dist=Normal, verbose=False)

grid_search = GridSearchCV(ngb, param_grid=param_grid, cv=5)
grid_search.fit(X_reg_train, Y_reg_train)
print(grid_search.best_params_)

/usr/local/lib/python3.7/site-packages/sklearn/model_selection/_search.py:814: DeprecationWarning: The default of the `iid` parameter will change from True to False in version 0.22 and will be removed in 0.24. This will change numeric results when test-set sizes are unequal.
  DeprecationWarning)

{'Base': DecisionTreeRegressor(criterion='friedman_mse', max_depth=4, max_features=None,
                      max_leaf_nodes=None, min_impurity_decrease=0.0,
                      min_impurity_split=None, min_samples_leaf=1,
                      min_samples_split=2, min_weight_fraction_leaf=0.0,
                      presort=False, random_state=None, splitter='best'), 'minibatch_frac': 1.0}