添加分布

首要任务是为你新的分布编写类。分布类必须继承适当的分布类型（要么是RegressionDistn，要么是ClassificationDistn），并且必须实现fit()和sample()方法。与分布兼容的分数应存储在一个名为score的类属性中，参数数量应存储在一个名为n_params的类属性中。类还必须将（内部）分布参数存储在_params实例属性中。此外，回归分布必须实现一个mean()方法来支持点预测。

我们将以拉普拉斯分布为例。拉普拉斯分布的概率密度函数为 $\frac{1}{2b} e^{-\frac{|x-\mu|}{b}}$，用户可见的参数为 $\mu \in \mathbb{R}$ 和 $b > 0$，我们将这些参数称为 loc 和 scale，以符合 scipy.stats 实现。

在NGBoost中，所有参数必须在内部表示为$\mathbb R$，因此我们需要将$(\mu, b)$重新参数化为$(\mu, \log(b))$。后者是我们在初始化Laplace对象和实现评分时需要处理的参数。

from scipy.stats import laplace as dist
import numpy as np
from ngboost.distns.distn import RegressionDistn
from ngboost.scores import LogScore

class LaplaceLogScore(LogScore): # will implement this later
    pass

class Laplace(RegressionDistn):

    n_params = 2
    scores = [LaplaceLogScore] # will implement this later

    def __init__(self, params):
        # save the parameters
        self._params = params

        # create other objects that will be useful later
        self.loc = params[0]
        self.logscale = params[1]
        self.scale = np.exp(params[1]) # since params[1] is log(scale)
        self.dist = dist(loc=self.loc, scale=self.scale)

    def fit(Y):
        m, s = dist.fit(Y) # use scipy's implementation
        return np.array([m, np.log(s)])

    def sample(self, m):
        return np.array([self.dist.rvs() for i in range(m)])

    def __getattr__(self, name): # gives us access to Laplace.mean() required for RegressionDist.predict()
        if name in dir(self.dist):
            return getattr(self.dist, name)
        return None

    @property
    def params(self):
        return {'loc':self.loc, 'scale':self.scale}

fit() 方法是一个类方法，它接受一个观测向量并拟合一个边缘分布。同时，sample() 应该从 $P(Y|X=x)$ 返回 $m$ 个样本，每个样本都是一个长度为 len(Y) 的向量。

这里我们利用了scipy.stats已经实现了拉普拉斯分布的事实，因此我们可以借用它的fit()方法，并在rvs()上添加一个薄包装来获取样本。我们还使用内部scipy.stats对象的__getattr__()方法来访问其mean()方法。

最后，我们编写了一个便捷方法 params()，当调用时，它会返回用户期望看到的分布参数，即 $(\mu, b)$，而不是 $(\mu, \log b)$。

为我们的分布实现评分

现在我们将注意力转向实现一个可以与此分布一起使用的评分。我们将以对数评分作为示例。

所有实现的分数应该继承适当的分数类并实现三个方法：

score() : 给定数据 Y 时，当前参数下的分数值
d_score() : 给定数据 Y 时，当前参数下得分的导数
metric() : 当前参数下的黎曼度量值

class LaplaceLogScore(LogScore):

    def score(self, Y):
        return -self.dist.logpdf(Y)

    def d_score(self, Y):
        D = np.zeros((len(Y), 2)) # first col is dS/d𝜇, second col is dS/d(log(b))
        D[:, 0] = np.sign(self.loc - Y)/self.scale
        D[:, 1] = 1 - np.abs(self.loc - Y)/self.scale
        return D

请注意，即使我们没有明确说明这些将是LaplaceLogScore类的属性，Laplace实例的属性也是使用self.attr符号引用的。当用户要求NGBoost使用LogScore的Laplace分布时，NGBoost将首先找到与Laplace兼容的对数分数实现，即LaplaceLogScore，并动态创建一个新类，该类既具有分布的属性，又具有适当的分数实现。为了使此功能正常工作，分布类Laplace必须具有一个scores类属性，其中包括实现LaplaceLogScore，并且LaplaceLogScore必须是LogScore的子类。只要满足这些条件，NGBoost就可以处理其余部分。

关于$\log b$和$\mu$的导数可以很容易地使用例如WolframAlpha来推导。

在这个例子中，我们不会费心去实现metric()，它会返回当前的费舍尔信息。原因是LogScore的NGBoost实现有一个默认的metric()方法，该方法使用蒙特卡罗方法通过gradient()方法和分布的sample()方法来近似费舍尔信息（这就是为什么我们需要实现sample()）。通过继承LogScore()，NGBoost不仅可以找到我们对拉普拉斯分布的实现，还可以回退到默认的metric()方法。更多内容稍后讨论。

将所有内容整合在一起：

class LaplaceLogScore(LogScore):

    def score(self, Y):
        return -self.dist.logpdf(Y)

    def d_score(self, Y):
        D = np.zeros((len(Y), 2)) # first col is dS/d𝜇, second col is dS/d(log(b))
        D[:, 0] = np.sign(self.loc - Y)/self.scale
        D[:, 1] = 1 - np.abs(self.loc - Y)/self.scale
        return D

class Laplace(RegressionDistn):

    n_params = 2
    scores = [LaplaceLogScore]

    def __init__(self, params):
        # save the parameters
        self._params = params

        # create other objects that will be useful later
        self.loc = params[0]
        self.logscale = params[1]
        self.scale = np.exp(params[1]) # since params[1] is log(scale)
        self.dist = dist(loc=self.loc, scale=self.scale)

    def fit(Y):
        m, s = dist.fit(Y) # use scipy's implementation
        return np.array([m, np.log(s)])

    def sample(self, m):
        return np.array([self.dist.rvs() for i in range(m)])

    def __getattr__(self, name): # gives us access to Laplace.mean() required for RegressionDist.predict()
        if name in dir(self.dist):
            return getattr(self.dist, name)
        return None

    @property
    def params(self):
        return {'loc':self.loc, 'scale':self.scale}

我们可以测试我们的方法：

from ngboost import NGBRegressor
from sklearn.datasets import load_boston
from sklearn.model_selection import train_test_split
from sklearn.metrics import mean_squared_error

X, Y = load_boston(True)
X_reg_train, X_reg_test, Y_reg_train, Y_reg_test = train_test_split(X, Y, test_size=0.2)

ngb = NGBRegressor(Dist=Laplace, Score=LogScore).fit(X_reg_train, Y_reg_train)
Y_preds = ngb.predict(X_reg_test)
Y_dists = ngb.pred_dist(X_reg_test)

# test Mean Squared Error
test_MSE = mean_squared_error(Y_preds, Y_reg_test)
print('Test MSE', test_MSE)

# test Negative Log Likelihood
test_NLL = -Y_dists.logpdf(Y_reg_test).mean()
print('Test NLL', test_NLL)

[iter 0] loss=3.5989 val_loss=0.0000 scale=1.0000 norm=6.8968
[iter 100] loss=2.8432 val_loss=0.0000 scale=1.0000 norm=5.3980
[iter 200] loss=2.4432 val_loss=0.0000 scale=1.0000 norm=3.5194
[iter 300] loss=2.2069 val_loss=0.0000 scale=1.0000 norm=2.4722
[iter 400] loss=2.0325 val_loss=0.0000 scale=1.0000 norm=1.9901
Test MSE 13.593387122611777
Test NLL 2.552651522788599

深入研究ngboost.distns的源代码以找到更多示例。如果您编写并测试了自己的分布，请通过提交拉取请求将其贡献给NGBoost！

审查分数

如果您实现了评分的审查版本，您可以使您的分布适合用于生存分析。score()、d_score() 和 metric() 方法的签名应该相同，但它们应该期望 Y 可以索引到两个数组，如 E, T = Y["Event"], Y["Time"]。此外，任何审查评分应通过一个名为 censored_scores 的类属性链接到分布类定义，而不是 scores。

由于截尾评分比其标准对应物更通用（完全观测数据是截尾数据的一个特例），如果你在NGBoost中实现了截尾评分，它将自动成为标准回归分析中可用的评分。无需单独实现回归评分或在scores类属性中注册它。

指标

正如我们所看到的，使用对数评分，作为开发者最简单的事情就是依赖默认的ngboost方法，该方法计算对数评分指标。

然而，与分布无关的默认方法速度较慢，因为它必须从分布中多次采样以构建度量的近似值。如果你想加快速度，那么你必须推导并实现特定于分布的黎曼度量，对于对数分数来说，这就是该分布的费舍尔信息矩阵。你必须根据内部ngboost参数化推导费舍尔（如果这与用户面对的参数化不同，例如$\log(\sigma)$，而不是$\sigma$）。推导费舍尔并不一定容易，因为你必须解析地计算期望值，但网上有很多推导费舍尔矩阵的例子可以参考。

例如，考虑Student's-t 分布。这个分布由自由度 $\nu$、均值 $\mu$ 和标准差 $\sigma$ 参数化。

这个分布的Fisher信息可以在这里找到，并且是

$$\begin{align*} \begin{bmatrix} \frac{\nu + 1}{(\nu + 3) \sigma^2} & 0 \\ 0 & \frac{\nu}{2(\nu + 3) \sigma^4} \end{bmatrix} \end{align*}$$

由于 $\sigma > 0$，对于NGBoost，我们必须用 $\log ( \sigma )$ 替换它。这需要我们重新参数化分布。要找到在这种重新参数化下的Fisher信息，我们可以按照Wikipedia上列出的步骤进行。

$\eta = (\mu, \sigma), \theta = (\mu, \log \sigma)$

$I_{\eta}(\eta) = J^T I_{\theta}(\theta) J$

其中 $J$ 是由定义的 $2 \times 2$ 雅可比矩阵

$$[J]_{ij} = \dfrac{\partial \theta_i}{\partial \eta_j}$$

其评估结果为

$$\begin{align*} J = J^T &= \begin{bmatrix} 1 & 0 \\ 0 & \frac{1}{\sigma} \end{bmatrix}\\ J^{-1} &= \begin{bmatrix} 1 & 0 \\ 0 & \sigma \end{bmatrix}\\ \end{align*}$$

因此，我们可以通过如下重新排列来获得所需的Fisher信息，

$$\begin{align*} I_{\theta}(\theta) &= J^{-1} I_{\eta}(\eta) J^{-1}\\ &= \begin{bmatrix} 1 & 0 \\ 0 & \sigma \end{bmatrix} \begin{bmatrix} \frac{\nu + 1}{(\nu + 3) \sigma^2} & 0 \\ 0 & \frac{\nu}{2(\nu + 3) \sigma^4} \end{bmatrix} \begin{bmatrix} 1 & 0 \\ 0 & \sigma \end{bmatrix}\\ &= \begin{bmatrix} \frac{\nu + 1}{(\nu + 3) \sigma^2} & 0 \\ 0 & \frac{\nu}{2(\nu + 3) \sigma^2} \end{bmatrix} \end{align*} $$

如果你不想使用对数评分（例如你想要CRP评分），那么ngboost（目前？）没有默认的方法来计算这个指标，你必须自己推导并实现它。这比推导Fisher信息更难，因为没有很多现成的例子。最一般的推导过程应该遵循这里的概述，用你想要的评分规则所诱导的散度（例如CRPS的L2）替换KL散度（由对数评分诱导），再次注意要针对内部的ngboost参数化进行推导，而不是面向用户的参数化。对于任何特定的评分，可能有一个特定的闭式表达式可以用来计算跨分布的指标（Fisher信息的表达式就是对数评分的用途），也可能没有——我实际上不知道这个问题的答案！但如果有的话，那可能暗示了该评分的metric()方法的某种默认实现。

添加分数

我们已经了解了如何为新的分布实现现有的评分，但在NGBoost中完全创建一个新的评分也很容易：只需创建一个继承Score的新类：

from ngboost.scores import Score

class SphericalScore(Score):
    pass

就是这样。这个分数的特定分布实现（例如LaplaceSphericalScore）应该继承SphericalScore。LogScore和CRPScore的实现可以在ngboost.scores中参考。