1.1. 线性模型

以下是一组用于回归的方法，其中目标值被期望为特征的线性组合。用数学符号表示，预测值 \(\hat{y}\) 可以写成

\[\hat{y}(w, x) = w_0 + w_1 x_1 + ... + w_p x_p\]

在整个模块中，我们将向量 \(w = (w_1, ..., w_p)\) 指定为 coef_，将 \(w_0\) 指定为 intercept_。

要使用广义线性模型进行分类，请参阅 逻辑回归。

1.1.1. 普通最小二乘法

LinearRegression 拟合一个系数为 \(w = (w_1, ..., w_p)\) 的线性模型，以最小化数据集中观测目标与线性近似预测目标之间的残差平方和。在数学上，它解决以下形式的问题：

\[\min_{w} || X w - y||_2^2\]

LinearRegression 在其 fit 方法中接收参数 X, y, sample_weight，并将线性模型的系数 \(w\) 存储在其 coef_ 和 intercept_ 属性中。

>>> from sklearn import linear_model
>>> reg = linear_model.LinearRegression()
>>> reg.fit([[0, 0], [1, 1], [2, 2]], [0, 1, 2])
LinearRegression()
>>> reg.coef_
array([0.5, 0.5])
>>> reg.intercept_
0.0

普通最小二乘法的系数估计依赖于特征的独立性。当特征相关且设计矩阵 \(X\) 的某些列具有近似的线性依赖关系时，设计矩阵变得接近奇异，结果是最小二乘估计对观测目标中的随机误差高度敏感，从而产生较大的方差。例如，当在没有实验设计的情况下收集数据时，可能会出现这种多重共线性情况。

示例

• 普通最小二乘法和岭回归

1.1.1.1. 非负最小二乘法

可以将所有系数限制为非负，这在它们代表某些物理或自然非负量（例如频率计数或商品价格）时非常有用。LinearRegression 接受一个布尔值 positive 参数：当设置为 True 时，应用 非负最小二乘法。

示例

• 非负最小二乘法

1.1.1.2. 普通最小二乘法的复杂度

最小二乘解是使用 \(X\) 的奇异值分解计算得出的。如果 \(X\) 是一个形状为 (n_samples, n_features) 的矩阵，且假设 \(n_{\text{samples}} \geq n_{\text{features}}\)，则该方法的计算成本为 \(O(n_{\text{samples}} n_{\text{features}}^2)\)。

1.1.2. 岭回归与分类

1.1.2.1. 回归

Ridge 回归通过对系数的大小施加惩罚来解决 普通最小二乘法 的一些问题。岭回归系数最小化惩罚后的残差平方和

\[\min_{w} || X w - y||_2^2 + \alpha ||w||_2^2\]

复杂度参数 \(\alpha \geq 0\) 控制收缩量：\(\alpha\) 的值越大，收缩量越大，因此系数对共线性的鲁棒性越强。

与其他线性模型一样，Ridge 在其 fit 方法中接收数组 X, y，并将线性模型的系数 \(w\) 存储在其 coef_ 成员中。

>>> from sklearn import linear_model
>>> reg = linear_model.Ridge(alpha=.5)
>>> reg.fit([[0, 0], [0, 0], [1, 1]], [0, .1, 1])
Ridge(alpha=0.5)
>>> reg.coef_
array([0.34545455, 0.34545455])
>>> reg.intercept_
np.float64(0.13636)

请注意，Ridge 类允许用户通过设置 solver="auto" 来自动选择求解器。当指定此选项时，Ridge 将在 "lbfgs", "cholesky" 和 "sparse_cg" 求解器之间进行选择。Ridge 将从上到下检查下表所示的条件。如果条件为真，则选择相应的求解器。

求解器	条件
‘lbfgs’	指定了 positive=True 选项。
‘cholesky’	输入数组 X 不是稀疏的。
‘sparse_cg’	以上条件均不满足。

示例

• 普通最小二乘法和岭回归

• 绘制岭系数作为正则化函数的函数

• 线性模型系数解释中的常见陷阱

• 岭系数作为 L2 正则化的函数

1.1.2.2. 分类

Ridge 回归器有一个分类器变体：RidgeClassifier。该分类器首先将二元目标转换为 {-1, 1}，然后将该问题视为回归任务，优化与上述相同的目标。预测类别对应于回归器预测值的符号。对于多类分类，问题被视为多输出回归，预测类别对应于具有最高值的输出。

使用（惩罚）最小二乘损失来拟合分类模型而不是更传统的逻辑损失或铰链损失，这似乎令人怀疑。然而，在实践中，所有这些模型在准确率或精确率/召回率方面都可能导致相似的交叉验证得分，而 RidgeClassifier 使用的惩罚最小二乘损失允许选择完全不同的数值求解器，这些求解器具有截然不同的计算性能曲线。

RidgeClassifier 可能比例如 LogisticRegression 在具有大量类别时快得多，因为它只需计算一次投影矩阵 \((X^T X)^{-1} X^T\)。

该分类器有时被称为具有线性核的 最小二乘支持向量机。

示例

• 使用稀疏特征对文本文档进行分类

1.1.2.3. 岭回归复杂度

该方法具有与 普通最小二乘法 相同阶数的复杂度。

1.1.2.4. 设置正则化参数：留一法交叉验证

RidgeCV 和 RidgeClassifierCV 实现了具有内置 alpha 参数交叉验证的岭回归/分类。它们的工作方式与 GridSearchCV 相同，只是默认使用高效的留一法 交叉验证。在使用默认 交叉验证 时，由于用于计算留一法误差的公式，alpha 不能为 0。详细信息请参阅 [RL2007]。

用法示例

>>> import numpy as np
>>> from sklearn import linear_model
>>> reg = linear_model.RidgeCV(alphas=np.logspace(-6, 6, 13))
>>> reg.fit([[0, 0], [0, 0.1], [1, 1]], [0, -0.1, 1])
RidgeCV(alphas=array([1.e-06, 1.e-05, 1.e-04, 1.e-03, 1.e-02, 1.e-01, 1.e+00, 1.e+01,
      1.e+02, 1.e+03, 1.e+04, 1.e+05, 1.e+06]))
>>> reg.alpha_
np.float64(0.1)

指定 cv 属性的值将触发使用 GridSearchCV 进行交叉验证，例如 cv=10 表示 10 折交叉验证，而不是留一法交叉验证。

参考文献

[RL2007]

“正则化最小二乘法笔记”，Rifkin & Lippert (技术报告, 课程幻灯片)。

1.1.3. Lasso

Lasso 是一种估计稀疏系数的线性模型，即它能够将系数精确设置为零。由于它倾向于偏好非零系数较少的解，从而有效地减少了给定解所依赖的特征数量，因此在某些上下文中非常有用。出于这个原因，Lasso 及其变体是压缩感知领域的基础。在某些条件下，它可以恢复精确的非零系数集（参见 压缩感知：基于 L1 先验的层析成像重构 (Lasso)）。

在数学上，它由一个带有额外正则化项的线性模型组成。需要最小化的目标函数为

\[\min_{w} P(w) = {\frac{1}{2n_{\text{samples}}} ||X w - y||_2 ^ 2 + \alpha ||w||_1}\]

因此，Lasso 估计量解决了添加了惩罚项 \(\alpha ||w||_1\) 的最小二乘法问题，其中 \(\alpha\) 是一个常数，\(||w||_1\) 是系数向量的 \(\ell_1\) 范数。

Lasso 类中的实现使用坐标下降算法来拟合系数。有关另一种实现，请参阅 最小角回归

>>> from sklearn import linear_model
>>> reg = linear_model.Lasso(alpha=0.1)
>>> reg.fit([[0, 0], [1, 1]], [0, 1])
Lasso(alpha=0.1)
>>> reg.predict([[1, 1]])
array([0.8])

lasso_path 函数对于较低级别的任务非常有用，因为它计算所有可能 alpha 值路径上的系数。

示例

• 用于稀疏信号的 L1 模型

• 压缩感知：使用 L1 先验（Lasso）进行断层扫描重建

• 线性模型系数解释中的常见陷阱

• Lasso 模型选择：AIC-BIC / 交叉验证

注意

使用 Lasso 进行特征选择

由于 Lasso 回归会产生稀疏模型，因此它可用于执行特征选择，详细信息请参见 基于 L1 的特征选择。

参考文献

以下参考文献解释了 Lasso 的起源以及 Lasso 问题的属性和用于收敛控制的对偶间隙计算。

• Robert Tibshirani. (1996) 通过 Lasso 进行回归收缩和选择。J. R. Stat. Soc. Ser. B Stat. Methodol., 58(1):267-288

• “大规模 L1 正则化最小二乘法的内点法”，S. J. Kim, K. Koh, M. Lustig, S. Boyd 和 D. Gorinevsky，发表于 IEEE 信号处理精选课题期刊，2007 (论文)

1.1.3.1. 带间隙安全筛选规则的坐标下降法

坐标下降法 (CD) 是一种求解最小化问题的策略，它每次只考虑一个特征 \(j\)。通过这种方式，优化问题被简化为更容易求解的一维问题

\[\min_{w_j} {\frac{1}{2n_{\text{samples}}} ||x_j w_j + X_{-j}w_{-j} - y||_2 ^ 2 + \alpha |w_j|}\]

索引 \(-j\) 表示除 \(j\) 之外的所有特征。解为

\[w_j = \frac{S(x_j^T (y - X_{-j}w_{-j}), \alpha)}{||x_j||_2^2}\]

其中软阈值函数 \(S(z, \alpha) = \operatorname{sign}(z) \max(0, |z|-\alpha)\)。注意，只要 \(\alpha \geq |z|\)，软阈值函数就完全为零。CD 求解器然后以循环方式遍历特征，按照 X 给出的顺序（selection="cyclic"）逐个选取特征，或者随机选取特征（selection="random"）。如果对偶间隙小于提供的公差 tol，它将停止。

数学细节

对偶间隙 \(G(w, v)\) 是 Lasso 当前原始目标函数 \(P(w)\) 与其最小值 \(P(w^\star)\) 之间差值的上界，即 \(G(w, v) \geq P(w) - P(w^\star)\)。它由 \(G(w, v) = P(w) - D(v)\) 给出，其中对偶目标函数为

\[D(v) = \frac{1}{2n_{\text{samples}}}(y^Tv - ||v||_2^2)\]

受约束于 \(v \in ||X^Tv||_{\infty} \leq n_{\text{samples}}\alpha\)。在最优状态下，对偶间隙为零，\(G(w^\star, v^\star) = 0\)（这被称为强对偶性）。使用（缩放后的）对偶变量 \(v = c r\)，当前残差 \(r = y - Xw\) 以及对偶缩放

\[\begin{split}c = \begin{cases} 1, & ||X^Tr||_{\infty} \leq n_{\text{samples}}\alpha, \\ \frac{n_{\text{samples}}\alpha}{||X^Tr||_{\infty}}, & \text{otherwise} \end{cases}\end{split}\]

停止准则为

\[\text{tol} \frac{||y||_2^2}{n_{\text{samples}}} < G(w, cr)\,.\]

加速坐标下降算法的一个巧妙方法是筛选特征，使得在最优解处 \(w_j = 0\)。间隙安全筛选规则就是这样一种工具。在优化算法执行过程中的任何位置，它们都可以告诉我们哪些特征可以安全地排除，即确定性地设为零。

参考文献

第一个参考文献解释了 scikit-learn 中使用的坐标下降求解器，其他参考文献则处理间隙安全筛选规则。

• Friedman, Hastie & Tibshirani. (2010). 通过坐标下降法进行广义线性模型的正则化路径。J Stat Softw 33(1), 1-22

• O. Fercoq, A. Gramfort, J. Salmon. (2015). 注意对偶间隙：Lasso 的更安全规则。机器学习研究论文集 37:333-342, 2015.

• E. Ndiaye, O. Fercoq, A. Gramfort, J. Salmon. (2017). 强制稀疏惩罚的间隙安全筛选规则。机器学习研究杂志 18(128):1-33, 2017.

1.1.3.2. 设置正则化参数

alpha 参数控制估计系数的稀疏程度。

1.1.3.2.1. 使用交叉验证

scikit-learn 提供了通过交叉验证设置 Lasso alpha 参数的对象：LassoCV 和 LassoLarsCV。LassoLarsCV 基于下面解释的 最小角回归 算法。

对于具有许多共线特征的高维数据集，LassoCV 通常更可取。然而，LassoLarsCV 的优点在于它可以探索更多相关的 alpha 参数值，并且如果样本数量相对于特征数量非常少，它通常比 LassoCV 快。

1.1.3.2.2. 基于信息准则的模型选择

或者，估计量 LassoLarsIC 建议使用赤池信息准则 (AIC) 和贝叶斯信息准则 (BIC)。这是一种在寻找最优 alpha 值时计算成本更低的选择，因为正则化路径只计算一次，而不是像使用 k 折交叉验证那样计算 k+1 次。

事实上，这些准则是在样本内训练集上计算的。简而言之，它们通过灵活性惩罚了不同 Lasso 模型过于乐观的得分（参考下文的“数学细节”部分）。

然而，此类准则需要对解的自由度进行适当估计，它们是针对大样本推导的（渐近结果），并假设所研究的候选模型是正确的模型。当问题条件不佳（例如特征多于样本）时，它们也往往会失效。

示例

• Lasso 模型选择：AIC-BIC / 交叉验证

• 通过信息准则进行 Lasso 模型选择

1.1.3.2.3. AIC 和 BIC 准则

AIC（以及 BIC）的定义在文献中可能有所不同。在本节中，我们提供了关于 scikit-learn 中计算准则的更多信息。

数学细节

AIC 准则定义为

\[AIC = -2 \log(\hat{L}) + 2 d\]

其中 \(\hat{L}\) 是模型的最大似然估计，\(d\) 是参数数量（在上一节中也称为自由度）。

BIC 的定义将常数 \(2\) 替换为 \(\log(N)\)

\[BIC = -2 \log(\hat{L}) + \log(N) d\]

其中 \(N\) 是样本数量。

对于线性高斯模型，最大对数似然定义为

\[\log(\hat{L}) = - \frac{n}{2} \log(2 \pi) - \frac{n}{2} \log(\sigma^2) - \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{2\sigma^2}\]

其中 \(\sigma^2\) 是噪声方差的估计值，\(y_i\) 和 \(\hat{y}_i\) 分别是真实目标值和预测目标值，\(n\) 是样本数量。

将最大对数似然代入 AIC 公式得到

\[AIC = n \log(2 \pi \sigma^2) + \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{\sigma^2} + 2 d\]

上述表达式的第一项有时会被舍弃，因为当提供 \(\sigma^2\) 时它是常数。此外，有时有人指出 AIC 等同于 \(C_p\) 统计量 [12]。然而，在严格意义上，它仅在差一个常数和乘法因子的前提下等同。

最后，我们上面提到 \(\sigma^2\) 是噪声方差的估计值。在 LassoLarsIC 中，当未提供 noise_variance 参数（默认情况）时，噪声方差通过无偏估计量进行估计 [13]，定义为

\[\sigma^2 = \frac{\sum_{i=1}^{n} (y_i - \hat{y}_i)^2}{n - p}\]

其中 \(p\) 是特征数量，\(\hat{y}_i\) 是使用普通最小二乘回归预测的目标值。注意，此公式仅在 n_samples > n_features 时有效。

References

[12]

Zou, Hui, Trevor Hastie, and Robert Tibshirani. “On the degrees of freedom of the lasso.” The Annals of Statistics 35.5 (2007): 2173-2192.

[13]

Cherkassky, Vladimir, and Yunqian Ma. “回归模型选择的比较。” 神经计算 15.7 (2003): 1691-1714。

1.1.3.2.4. 与 SVM 正则化参数的比较

alpha 与 SVM 正则化参数 C 之间的等价关系由 alpha = 1 / C 或 alpha = 1 / (n_samples * C) 给出，具体取决于估计器和模型优化的确切目标函数。

1.1.4. 多任务 Lasso

MultiTaskLasso 是一种线性模型，它联合估计多个回归问题的稀疏系数：y 是一个形状为 (n_samples, n_tasks) 的二维数组。其约束条件是所选特征在所有回归问题（也称为任务）中是相同的。

下图比较了通过简单 Lasso 或 MultiTaskLasso 获得的系数矩阵 W 中非零条目的位置。Lasso 估计产生零星的非零值，而 MultiTaskLasso 的非零值是完整的列。

拟合时间序列模型，要求任何活跃特征在所有时间点都保持活跃。

示例

• 使用多任务 Lasso 进行联合特征选择

数学细节

在数学上，它由一个使用混合 \(\ell_1\) \(\ell_2\)-范数进行正则化的线性模型组成。需要最小化的目标函数为

\[\min_{W} { \frac{1}{2n_{\text{samples}}} ||X W - Y||_{\text{Fro}} ^ 2 + \alpha ||W||_{21}}\]

其中 \(\text{Fro}\) 表示 Frobenius 范数

\[||A||_{\text{Fro}} = \sqrt{\sum_{ij} a_{ij}^2}\]

且 \(\ell_1\) \(\ell_2\) 读取为

\[||A||_{2 1} = \sum_i \sqrt{\sum_j a_{ij}^2}.\]

MultiTaskLasso 类中的实现使用坐标下降算法来拟合系数。

1.1.5. Elastic-Net

ElasticNet 是一种线性回归模型，使用系数的 \(\ell_1\) 和 \(\ell_2\)-范数正则化进行训练。这种组合允许学习一个像 Lasso 那样权重稀疏（非零值较少）的模型，同时仍保持 Ridge 的正则化属性。我们使用 l1_ratio 参数控制 \(\ell_1\) 和 \(\ell_2\) 的凸组合。

当存在多个彼此相关的特征时，Elastic-net 很有用。Lasso 可能会随机选择其中一个，而 elastic-net 更有可能同时选择它们。

在 Lasso 和 Ridge 之间进行权衡的一个实际优点是，它允许 Elastic-Net 继承一些 Ridge 在旋转下的稳定性。

在这种情况下，需要最小化的目标函数为

\[\min_{w} { \frac{1}{2n_{\text{samples}}} ||X w - y||_2 ^ 2 + \alpha \rho ||w||_1 + \frac{\alpha(1-\rho)}{2} ||w||_2 ^ 2}\]

ElasticNetCV 类可用于通过交叉验证设置参数 alpha (\(\alpha\)) 和 l1_ratio (\(\rho\))。

示例

• 用于稀疏信号的 L1 模型

• Lasso、Lasso-LARS 和 Elastic Net 路径

• 使用预计算 Gram 矩阵和加权样本拟合 Elastic Net

参考文献

以下两个参考文献解释了 scikit-learn 坐标下降求解器中使用的迭代，以及用于收敛控制的对偶间隙计算。

• “通过坐标下降法进行广义线性模型的正则化路径”，Friedman, Hastie & Tibshirani, J Stat Softw, 2010 (论文)。

• “大规模 L1 正则化最小二乘法的内点法”，S. J. Kim, K. Koh, M. Lustig, S. Boyd 和 D. Gorinevsky，发表于 IEEE 信号处理精选课题期刊，2007 (论文)

1.1.6. 多任务 Elastic-Net

MultiTaskElasticNet 是一种 elastic-net 模型，它联合估计多个回归问题的稀疏系数：Y 是形状为 (n_samples, n_tasks) 的 2D 数组。其约束条件是所选特征在所有回归问题（也称为任务）中是相同的。

在数学上，它由一个使用混合 \(\ell_1\) \(\ell_2\)-范数和 \(\ell_2\)-范数进行正则化的线性模型组成。需要最小化的目标函数为

\[\min_{W} { \frac{1}{2n_{\text{samples}}} ||X W - Y||_{\text{Fro}}^2 + \alpha \rho ||W||_{2 1} + \frac{\alpha(1-\rho)}{2} ||W||_{\text{Fro}}^2}\]

MultiTaskElasticNet 类中的实现使用坐标下降算法来拟合系数。

MultiTaskElasticNetCV 类可用于通过交叉验证设置参数 alpha (\(\alpha\)) 和 l1_ratio (\(\rho\))。

1.1.7. 最小角回归

最小角回归 (LARS) 是一种用于高维数据的回归算法，由 Bradley Efron, Trevor Hastie, Iain Johnstone 和 Robert Tibshirani 开发。LARS 类似于前向逐步回归。在每一步，它都会找到与目标相关性最高的特征。当存在多个具有相同相关性的特征时，它不会继续沿同一特征前进，而是沿这些特征之间的等角方向前进。

LARS 的优点是

• 在特征数量显著大于样本数量的上下文中，它在数值上是有效的。

• 在计算上，它与前向选择一样快，并且与普通最小二乘法具有相同阶数的复杂度。

• 它产生完整的逐段线性解路径，这在交叉验证或类似调整模型的尝试中很有用。

• 如果两个特征与目标的相关性几乎相同，那么它们的系数应该以大致相同的速率增加。因此，该算法表现得符合直觉，而且更稳定。

• 它可以很容易地修改以产生其他估计器的解，例如 Lasso。

LARS 方法的缺点包括

• 由于 LARS 基于残差的迭代重拟合，它似乎对噪声的影响特别敏感。Weisberg 在 Efron 等人 (2004) Annals of Statistics 文章的讨论部分详细讨论了这个问题。

LARS 模型可以通过估计器 Lars，或者其低级实现 lars_path 或 lars_path_gram 来使用。

1.1.8. LARS Lasso

LassoLars 是一种使用 LARS 算法实现的 lasso 模型，与基于坐标下降的实现不同，这产生的是精确解，它作为其系数范数的函数是逐段线性的。

>>> from sklearn import linear_model
>>> reg = linear_model.LassoLars(alpha=.1)
>>> reg.fit([[0, 0], [1, 1]], [0, 1])
LassoLars(alpha=0.1)
>>> reg.coef_
array([0.6, 0.        ])

示例

• Lasso、Lasso-LARS 和 Elastic Net 路径

LARS 算法几乎免费提供沿正则化参数的系数完整路径，因此一种常见操作是使用函数 lars_path 或 lars_path_gram 检索该路径。

数学公式

该算法类似于前向逐步回归，但不是在每一步包含特征，而是估计系数沿每个特征与残差相关性的等角方向增加。

LARS 解不是给出一个向量结果，而是由一条曲线组成，表示参数向量 \(\ell_1\) 范数的每个值的解。完整系数路径存储在形状为 (n_features, max_features + 1) 的数组 coef_path_ 中。第一列始终为零。

References

• 原始算法在 Hastie 等人的论文 最小角回归 中有详细说明。

1.1.9. 正交匹配追踪 (OMP)

OrthogonalMatchingPursuit 和 orthogonal_mp 实现了 OMP 算法，用于逼近对非零系数数量（即 \(\ell_0\) 伪范数）施加约束的线性模型拟合。

作为一种类似于 最小角回归 的前向特征选择方法，正交匹配追踪可以用固定数量的非零元素逼近最优解向量

\[\underset{w}{\operatorname{arg\,min\,}} ||y - Xw||_2^2 \text{ subject to } ||w||_0 \leq n_{\text{nonzero_coefs}}\]

或者，正交匹配追踪可以针对特定误差而不是特定数量的非零系数。这可以表示为

\[\underset{w}{\operatorname{arg\,min\,}} ||w||_0 \text{ subject to } ||y-Xw||_2^2 \leq \text{tol}\]

OMP 基于一种贪婪算法，该算法在每一步都包含与当前残差相关性最高的原子。它类似于更简单的匹配追踪 (MP) 方法，但更好之处在于，在每次迭代中，残差都是使用先前选择的字典元素空间上的正交投影重新计算的。

示例

• 正交匹配追踪

参考文献

• https://www.cs.technion.ac.il/~ronrubin/Publications/KSVD-OMP-v2.pdf

• 具有时频字典的匹配追踪, S. G. Mallat, Z. Zhang, 1993.

1.1.10. 贝叶斯回归

贝叶斯回归技术可用于将正则化参数包含在估计过程中：正则化参数不是硬性设定的，而是根据手头的数据进行调整。

这可以通过在模型超参数上引入 无信息先验 来实现。岭回归与分类 中使用的 \(\ell_{2}\) 正则化等同于在系数 \(w\) 上使用精度为 \(\lambda^{-1}\) 的高斯先验下找到最大后验估计。与其手动设置 lambda，不如将其视为要从数据中估计的随机变量。

为了获得一个完全概率化的模型，假设输出 \(y\) 是围绕 \(X w\) 的高斯分布

\[p(y|X,w,\alpha) = \mathcal{N}(y|X w,\alpha^{-1})\]

其中 \(\alpha\) 再次被视为要从数据中估计的随机变量。

贝叶斯回归的优点是

• 它适应手头的数据。

• 它可用于将正则化参数包含在估计过程中。

贝叶斯回归的缺点包括

• 模型的推断可能非常耗时。

参考文献

• C. Bishop: 模式识别与机器学习 中给出了贝叶斯方法的良好入门。

• 原始算法在 Radford M. Neal 所著的 神经网络的贝叶斯学习 一书中进行了详细说明。

1.1.10.1. 贝叶斯岭回归

BayesianRidge 估计上述回归问题的概率模型。系数 \(w\) 的先验由球形高斯分布给出

\[p(w|\lambda) = \mathcal{N}(w|0,\lambda^{-1}\mathbf{I}_{p})\]

\(\alpha\) 和 \(\lambda\) 的先验被选择为 伽马分布，这是高斯精度共轭先验。生成的模型被称为贝叶斯岭回归，类似于经典的 Ridge。

参数 \(w\), \(\alpha\) 和 \(\lambda\) 在模型拟合期间联合估计，正则化参数 \(\alpha\) 和 \(\lambda\) 通过最大化对数边际似然进行估计。scikit-learn 实现基于 (Tipping, 2001) 附录 A 中描述的算法，其中参数 \(\alpha\) 和 \(\lambda\) 的更新是按照 (MacKay, 1992) 的建议进行的。最大化过程的初始值可以使用超参数 alpha_init 和 lambda_init 设置。

还有四个超参数，\(\alpha_1\), \(\alpha_2\), \(\lambda_1\) 和 \(\lambda_2\)，它们是 \(\alpha\) 和 \(\lambda\) 上伽马先验分布的参数。这些通常被选择为无信息的。默认情况下 \(\alpha_1 = \alpha_2 = \lambda_1 = \lambda_2 = 10^{-6}\)。

贝叶斯岭回归用于回归

>>> from sklearn import linear_model
>>> X = [[0., 0.], [1., 1.], [2., 2.], [3., 3.]]
>>> Y = [0., 1., 2., 3.]
>>> reg = linear_model.BayesianRidge()
>>> reg.fit(X, Y)
BayesianRidge()

拟合后，模型可用于预测新值

>>> reg.predict([[1, 0.]])
array([0.50000013])

模型系数 \(w\) 可以被访问

>>> reg.coef_
array([0.49999993, 0.49999993])

由于贝叶斯框架，找到的权重与 普通最小二乘法 找到的权重略有不同。然而，贝叶斯岭回归对不适定问题更稳健。

示例

• 使用贝叶斯岭回归进行曲线拟合

参考文献

• Christopher M. Bishop: 模式识别与机器学习, 2006, 3.3 节

• David J. C. MacKay, 贝叶斯插值, 1992.

• Michael E. Tipping, 稀疏贝叶斯学习与相关向量机, 2001.

1.1.10.2. 自动相关性确定 - ARD

自动相关性确定（在 ARDRegression 中实现）是一种类似于 贝叶斯岭回归 的线性模型，但它会产生更稀疏的系数 \(w\) [1] [2]。

ARDRegression 在 \(w\) 上设置了不同的先验：它放弃了球形高斯分布，转而使用中心椭圆高斯分布。这意味着每个系数 \(w_{i}\) 本身可以从高斯分布中抽取，该分布以零为中心，精度为 \(\lambda_{i}\)

\[p(w|\lambda) = \mathcal{N}(w|0,A^{-1})\]

其中 \(A\) 是正定对角矩阵，且 \(\text{diag}(A) = \lambda = \{\lambda_{1},...,\lambda_{p}\}\)。

与 贝叶斯岭回归 相比，\(w_{i}\) 的每个坐标都有自己的标准差 \(\frac{1}{\lambda_i}\)。所有 \(\lambda_i\) 上的先验被选择为由超参数 \(\lambda_1\) 和 \(\lambda_2\) 给出的相同伽马分布。

ARD 在文献中也称为稀疏贝叶斯学习和相关向量机 [3] [4]。

有关 ARD 与 贝叶斯岭回归 之间的深入比较，请参阅 比较线性贝叶斯回归器。

有关 Lasso, ARD 和 ElasticNet 等各种方法在相关数据上的比较，请参阅 用于稀疏信号的 L1 模型。

References

[1]

Christopher M. Bishop: 模式识别与机器学习, 第 7.2.1 章

[2]

David Wipf and Srikantan Nagarajan: 自动相关性确定的新视角

[3]

Michael E. Tipping: 稀疏贝叶斯学习与相关向量机

[4]

Tristan Fletcher: 相关向量机解释

1.1.11. 逻辑回归

逻辑回归在 LogisticRegression 中实现。尽管名称如此，但在 scikit-learn/ML 命名法中，它是作为分类的线性模型而非回归模型实现的。逻辑回归在文献中也称为 logit 回归、最大熵分类 (MaxEnt) 或对数线性分类器。在该模型中，描述单次试验可能结果的概率是使用 逻辑函数 进行建模的。

此实现可以拟合二元、一对多 (One-vs-Rest) 或多项逻辑回归，并带有可选的 \(\ell_1\), \(\ell_2\) 或 Elastic-Net 正则化。

注意

正则化

默认应用正则化，这在机器学习中很常见，但在统计学中则不然。正则化的另一个优点是它提高了数值稳定性。不进行正则化等同于将 C 设置为一个非常大的值。

注意

逻辑回归作为广义线性模型 (GLM) 的特例

逻辑回归是 广义线性模型 的一个特例，具有二项式/伯努利条件分布和 Logit 链接。逻辑回归的数值输出（即预测概率）可以通过应用阈值（默认为 0.5）作为分类器使用。这就是它在 scikit-learn 中的实现方式，因此它期望分类目标，使逻辑回归成为分类器。

示例

• 逻辑回归中的 L1 惩罚和稀疏性

• L1 逻辑回归的正则化路径

• 多项式和一对多逻辑回归的决策边界

• 在 20newgroups 上进行多类稀疏逻辑回归

• 使用多项式逻辑 + L1 进行 MNIST 分类

• 绘制分类概率

1.1.11.1. 二元情况

为了方便起见，我们假设数据点 \(i\) 的目标 \(y_i\) 取集合 \(\{0, 1\}\) 中的值。一旦拟合，LogisticRegression 的 predict_proba 方法预测正类的概率 \(P(y_i=1|X_i)\) 为

\[\hat{p}(X_i) = \operatorname{expit}(X_i w + w_0) = \frac{1}{1 + \exp(-X_i w - w_0)}.\]

作为优化问题，具有正则化项 \(r(w)\) 的二元类逻辑回归最小化以下代价函数

(1)\[\min_{w} \frac{1}{S}\sum_{i=1}^n s_i \left(-y_i \log(\hat{p}(X_i)) - (1 - y_i) \log(1 - \hat{p}(X_i))\right) + \frac{r(w)}{S C}\,,\]

其中 \({s_i}\) 对应于用户为特定训练样本分配的权重（向量 \(s\) 由类别权重和样本权重的逐元素乘法形成），并且求和 \(S = \sum_{i=1}^n s_i\)。

我们目前通过参数 C 和 l1_ratio 为正则化或惩罚项 \(r(w)\) 提供四种选择

惩罚	\(r(w)\)
none (C=np.inf)	\(0\)
\(\ell_1\) (l1_ratio=1)	\(\|w\|_1\)
\(\ell_2\) (l1_ratio=0)	\(\frac{1}{2}\|w\|_2^2 = \frac{1}{2}w^T w\)
ElasticNet (0<l1_ratio<1)	\(\frac{1 - \rho}{2}w^T w + \rho \|w\|_1\)

对于 ElasticNet，\(\rho\)（对应于 l1_ratio 参数）控制 \(\ell_1\) 正则化与 \(\ell_2\) 正则化的强度。当 \(\rho = 1\) 时，Elastic-Net 等同于 \(\ell_1\)，当 \(\rho=0\) 时等同于 \(\ell_2\)。

请注意，类别权重和样本权重的比例将影响优化问题。例如，将样本权重乘以常数 \(b>0\) 等同于将（逆）正则化强度 C 乘以 \(b\)。

1.1.11.2. 多项情况

二元情况可以扩展到 \(K\) 个类别，从而得出多项逻辑回归，另请参阅 对数线性模型。

注意

可以通过利用所有类别概率必须相加为 1 的事实，仅使用 \(K-1\) 个权重向量对 \(K\) 类分类模型进行参数化。为了实现简便以及保持关于类别排序的对称归纳偏置，我们故意选择使用 \(K\) 个权重向量对模型进行过参数化，见 [16]。当使用正则化时，这种影响变得尤为重要。过参数化的选择对于未惩罚模型可能是有害的，因为正如 [16] 中所示，解可能不是唯一的。

数学细节

令 \(y_i \in \{1, \ldots, K\}\) 为观察值 \(i\) 的标签（序数）编码目标变量。我们现在拥有的不是单个系数向量，而是一个系数矩阵 \(W\)，其中每个行向量 \(W_k\) 对应于类别 \(k\)。我们旨在通过 predict_proba 预测类别概率 \(P(y_i=k|X_i)\) 为

\[\hat{p}_k(X_i) = \frac{\exp(X_i W_k + W_{0, k})}{\sum_{l=0}^{K-1} \exp(X_i W_l + W_{0, l})}.\]

优化的目标变为

\[\min_W -\frac{1}{S}\sum_{i=1}^n \sum_{k=0}^{K-1} s_{ik} [y_i = k] \log(\hat{p}_k(X_i)) + \frac{r(W)}{S C}\,,\]

其中 \([P]\) 表示 Iverson 括号，如果 \(P\) 为假则求值为 \(0\)，否则求值为 \(1\)。

同样，\(s_{ik}\) 是用户分配的权重（样本权重和类别权重的乘积），其总和为 \(S = \sum_{i=1}^n \sum_{k=0}^{K-1} s_{ik}\)。

我们目前通过参数 C 和 l1_ratio 为正则化或惩罚项 \(r(W)\) 提供四种选择，其中 \(m\) 是特征数量

惩罚	\(r(W)\)
none (C=np.inf)	\(0\)
\(\ell_1\) (l1_ratio=1)	\(\|W\|_{1,1} = \sum_{i=1}^m\sum_{j=1}^{K}|W_{i,j}|\)
\(\ell_2\) (l1_ratio=0)	\(\frac{1}{2}\|W\|_F^2 = \frac{1}{2}\sum_{i=1}^m\sum_{j=1}^{K} W_{i,j}^2\)
ElasticNet (0<l1_ratio<1)	\(\frac{1 - \rho}{2}\|W\|_F^2 + \rho \|W\|_{1,1}\)

1.1.11.3. 求解器

LogisticRegression 类中实现的求解器有 “lbfgs”, “liblinear”, “newton-cg”, “newton-cholesky”, “sag” 和 “saga”

下表总结了每个求解器支持的惩罚项和多项多类分类

	求解器
惩罚项	‘lbfgs’	‘liblinear’		‘newton-cg’	‘newton-cholesky’	‘sag’	‘saga’
L2 惩罚	是	是		是	是	是	是
L1 惩罚	否	是		否	否	否	是
Elastic-Net (L1 + L2)	否	否		否	否	否	是
无惩罚	是	否		是	是	是	是
多类支持
多项式多类	是	否		是	是	是	是
行为
惩罚截距（坏）	否	是		否	否	否	否
大型数据集更快	否	否		否	否	是	是
对未缩放数据集稳健	是	是		是	是	否	否

“lbfgs” 求解器因其稳健性而被默认使用。对于 n_samples >> n_features，"newton-cholesky" 是一个不错的选择，可以达到高精度（极小的 tol 值）。对于大型数据集，"saga" 求解器通常比 "lbfgs" 更快，尤其是在低精度（高 tol）的情况下。对于大型数据集，您还可以考虑使用带有 loss="log_loss" 的 SGDClassifier，它可能更快但需要更多的调整。

1.1.11.3.1. 求解器之间的差异

当 fit_intercept=False 且拟合的 coef_（或者）要预测的数据为零时，使用 solver=liblinear 的 LogisticRegression 或 LinearSVC 与直接使用外部 liblinear 库之间可能存在得分差异。这是因为对于 decision_function 为零的样本，LogisticRegression 和 LinearSVC 预测负类，而 liblinear 预测正类。请注意，具有 fit_intercept=False 且具有许多 decision_function 为零的样本的模型很可能是一个欠拟合的坏模型，建议您设置 fit_intercept=True 并增加 intercept_scaling。

求解器详情

• 求解器 "liblinear" 使用坐标下降 (CD) 算法，并依赖于随 scikit-learn 一起发布的出色的 C++ LIBLINEAR 库。然而，在 liblinear 中实现的 CD 算法无法学习真正的多项式（多类）模型。如果您仍然想在多类问题上使用 "liblinear"，可以使用“一对多”方案 OneVsRestClassifier(LogisticRegression(solver="liblinear"))，见 :class:`~sklearn.multiclass.OneVsRestClassifier。注意，与“一对多”方案相比，最小化多项式损失预计会给出校准更好的结果。对于 \(\ell_1\) 正则化，sklearn.svm.l1_min_c 允许计算 C 的下界，以获得非“空”（所有特征权重为零）的模型。

• “lbfgs”, “newton-cg”, “newton-cholesky” 和 “sag” 求解器仅支持 \(\ell_2\) 正则化或无正则化，并且被发现对于某些高维数据收敛得更快。这些求解器（以及 "saga"）学习真正的多项逻辑回归模型 [5]。

• "sag" 求解器使用随机平均梯度下降 [6]。当样本数量和特征数量都很大时，它对于大型数据集比其他求解器更快。

• "saga" 求解器 [7] 是 "sag" 的变体，也支持非平滑 \(\ell_1\) 惩罚（l1_ratio=1）。因此，这是稀疏多项逻辑回归的首选求解器。它也是唯一支持 Elastic-Net (0 < l1_ratio < 1) 的求解器。

• "lbfgs" 是一种优化算法，它近似 Broyden–Fletcher–Goldfarb–Shanno 算法 [8]，该算法属于拟牛顿法。因此，它可以处理各种不同的训练数据，因此是默认求解器。然而，它的性能在缩放不良的数据集以及具有罕见类别的独热编码分类特征的数据集上会受到影响。

• "newton-cholesky" 求解器是一种精确的牛顿求解器，它计算 Hessian 矩阵并求解生成的线性系统。对于 n_samples >> n_features 来说，这是一个非常好的选择，可以达到高精度（极小的 tol 值），但有一些缺点：仅支持 \(\ell_2\) 正则化。此外，由于显式计算了 Hessian 矩阵，内存使用与 n_features 以及 n_classes 呈二次方依赖。

有关其中一些求解器的比较，请参阅 [9]。

References

[5]

Christopher M. Bishop: 模式识别与机器学习, 第 4.3.4 章

[6]

Mark Schmidt, Nicolas Le Roux, and Francis Bach: 用随机平均梯度最小化有限和。

[7]

Aaron Defazio, Francis Bach, Simon Lacoste-Julien: SAGA：一种支持非强凸复合目标的快速增量梯度方法。

[8]

https://en.wikipedia.org/wiki/Broyden%E2%80%93Fletcher%E2%80%93Goldfarb%E2%80%93Shanno_algorithm

[9]

Thomas P. Minka “逻辑回归数值优化器的比较”

[16] (1,2)

Simon, Noah, J. Friedman and T. Hastie. “用于组惩罚多响应和多项回归的分块下降算法。”

注意

使用稀疏逻辑回归进行特征选择

带有 \(\ell_1\) 惩罚的逻辑回归会产生稀疏模型，因此可用于执行特征选择，详细信息请参见 基于 L1 的特征选择。

注意

P 值估计

在无惩罚回归的情况下，可以获得系数的 p 值和置信区间。statsmodels 包 原生支持这一点。在 sklearn 中，也可以改用自助法 (bootstrapping)。

LogisticRegressionCV 实现了具有内置交叉验证支持的逻辑回归，以便根据 scoring 属性找到最优的 C 和 l1_ratio 参数。由于热启动（参见 词汇表），发现 "newton-cg", "sag", "saga" 和 "lbfgs" 求解器对于高维稠密数据更快。

1.1.12. 广义线性模型

广义线性模型 (GLM) 在两个方面扩展了线性模型 [10]。首先，预测值 \(\hat{y}\) 通过反链接函数 \(h\) 与输入变量 \(X\) 的线性组合相关联：

\[\hat{y}(w, X) = h(Xw)。\]

其次，平方损失函数被指数族分布（更准确地说是再生指数离散模型 (EDM) [11]）的单位偏差 \(d\) 所取代。

最小化问题变为：

\[\min_{w} \frac{1}{2 n_{\text{samples}}} \sum_i d(y_i, \hat{y}_i) + \frac{\alpha}{2} ||w||_2^2,\]

其中 \(\alpha\) 是 L2 正则化惩罚项。当提供样本权重时，平均值将变为加权平均值。

下表列出了一些特定的 EDM 及其单位偏差：

分布	目标域	单位偏差 \(d(y, \hat{y})\)
正态分布	\(y \in (-\infty, \infty)\)	\((y-\hat{y})^2\)
伯努利分布	\(y \in \{0, 1\}\)	\(2({y}\log\frac{y}{\hat{y}}+({1}-{y})\log\frac{{1}-{y}}{{1}-\hat{y}})\)
分类分布	\(y \in \{0, 1, ..., k\}\)	\(2\sum_{i \in \{0, 1, ..., k\}} I(y = i) y_\text{i}\log\frac{I(y = i)}{\hat{I(y = i)}}\)
泊松分布	\(y \in [0, \infty)\)	\(2(y\log\frac{y}{\hat{y}}-y+\hat{y})\)
伽马分布	\(y \in (0, \infty)\)	\(2(\log\frac{\hat{y}}{y}+\frac{y}{\hat{y}}-1)\)
逆高斯分布	\(y \in (0, \infty)\)	\(\frac{(y-\hat{y})^2}{y\hat{y}^2}\)

这些分布的概率密度函数 (PDF) 如下图所示：

泊松分布、Tweedie 分布 (power=1.5) 和伽马分布在不同均值 (\(\mu\)) 下的随机变量 Y 的 PDF。观察到泊松分布和 Tweedie 分布 (power=1.5) 在 \(Y=0\) 处有点质量，而目标域严格为正的伽马分布则没有。

伯努利分布是一种离散概率分布，用于建模伯努利试验——即只有两种互斥结果的事件。分类分布是伯努利分布在分类随机变量上的推广。伯努利分布中的随机变量有两个可能的结果，而分类随机变量可以取 K 个可能类别中的一个，且每个类别的概率是分别指定的。

分布的选择取决于具体问题：

• 如果目标值 \(y\) 是计数（非负整数值）或相对频率（非负），您可以使用带有对数链接的泊松分布。

• 如果目标值是正值且偏斜，您可以尝试带有对数链接的伽马分布。

• 如果目标值的尾部比伽马分布更重，您可以尝试逆高斯分布（或 Tweedie 分布中更高方差的幂）。

• 如果目标值 \(y\) 是概率，您可以使用伯努利分布。带 logit 链接的伯努利分布可用于二分类。带 softmax 链接的分类分布可用于多分类。

使用案例示例...

• 农业/天气建模：每年降雨事件次数（泊松分布）、每次降雨量（伽马分布）、每年总降雨量（Tweedie / 复合泊松伽马分布）。

• 风险建模/保险费定价：每年每位投保人的索赔事件次数（泊松分布）、每次事件成本（伽马分布）、每年每位投保人的总成本（Tweedie / 复合泊松伽马分布）。

• 信用违约：贷款无法偿还的概率（伯努利分布）。

• 欺诈检测：金融交易（如现金转账）为欺诈交易的概率（伯努利分布）。

• 预测性维护：每年生产中断事件次数（泊松分布）、中断持续时间（伽马分布）、每年总中断时间（Tweedie / 复合泊松伽马分布）。

• 医疗药物测试：一系列试验中治愈患者的概率或患者产生副作用的概率（伯努利分布）。

• 新闻分类：将新闻文章分类为三个类别，即商业新闻、政治新闻和娱乐新闻（分类分布）。

References

[10]

McCullagh, Peter; Nelder, John (1989). Generalized Linear Models, Second Edition. Boca Raton: Chapman and Hall/CRC. ISBN 0-412-31760-5.

[11]

Jørgensen, B. (1992). The theory of exponential dispersion models and analysis of deviance. Monografias de matemática, no. 51. 另请参阅 Exponential dispersion model.

1.1.12.1. 用法

TweedieRegressor 实现了 Tweedie 分布的广义线性模型，通过适当的 power 参数，可以对上述任何分布进行建模。具体来说：

• power = 0：正态分布。在这种情况下，通常使用 Ridge 或 ElasticNet 等特定估计器更为合适。

• power = 1：泊松分布。为了方便起见，提供了 PoissonRegressor。它严格等同于 TweedieRegressor(power=1, link='log')。

• power = 2：伽马分布。为了方便起见，提供了 GammaRegressor。它严格等同于 TweedieRegressor(power=2, link='log')。

• power = 3：逆高斯分布。

链接函数由 link 参数确定。

用法示例

>>> from sklearn.linear_model import TweedieRegressor
>>> reg = TweedieRegressor(power=1, alpha=0.5, link='log')
>>> reg.fit([[0, 0], [0, 1], [2, 2]], [0, 1, 2])
TweedieRegressor(alpha=0.5, link='log', power=1)
>>> reg.coef_
array([0.2463, 0.4337])
>>> reg.intercept_
np.float64(-0.7638)

示例

• 泊松回归和非正态损失

• 保险索赔上的 Tweedie 回归

实际考量...

在拟合之前，应对特征矩阵 X 进行标准化。这确保了正则化项对所有特征一视同仁。

由于线性预测变量 \(Xw\) 可能为负，而泊松分布、伽马分布和逆高斯分布不支持负值，因此必须应用能保证非负性的反链接函数。例如，使用 link='log' 时，反链接函数变为 \(h(Xw)=\exp(Xw)\)。

如果您想对相对频率（即每个暴露量（时间、体积等）的计数）进行建模，可以使用泊松分布，并将 \(y=\frac{\mathrm{counts}}{\mathrm{exposure}}\) 作为目标值，并将 \(\mathrm{exposure}\) 作为样本权重。具体示例请参阅 保险索赔的 Tweedie 回归。

在对 TweedieRegressor 的 power 参数进行交叉验证时，建议指定显式的 scoring 函数，因为默认的评分器 TweedieRegressor.score 本身就是 power 的函数。

1.1.13. 随机梯度下降 - SGD

随机梯度下降是一种简单且非常高效的线性模型拟合方法。当样本数量（和特征数量）非常大时，它特别有用。partial_fit 方法允许在线学习/外存学习。

SGDClassifier 和 SGDRegressor 类提供了使用不同（凸）损失函数和不同正则化项来拟合分类和回归线性模型的功能。例如，使用 loss="log"，SGDClassifier 拟合逻辑回归模型，而使用 loss="hinge" 则拟合线性支持向量机 (SVM)。

您可以参考专门的 随机梯度下降 文档章节了解更多详细信息。

1.1.13.1. 感知机 (Perceptron)

Perceptron 是另一种适用于大规模学习的简单分类算法，它衍生自 SGD。默认情况下：

• 它不需要学习率。

• 它没有被正则化（惩罚）。

• 它仅在出错时更新模型。

最后一个特征意味着感知机的训练速度略快于使用 hinge 损失的 SGD，并且产生的模型更稀疏。

实际上，Perceptron 是 SGDClassifier 类的包装器，使用了感知机损失和恒定学习率。有关更多详细信息，请参阅 SGD 过程的 数学部分。

1.1.13.2. 被动攻击算法 (Passive Aggressive Algorithms)

被动攻击 (PA) 算法是另一种衍生自 SGD 的大规模在线学习算法系列（PA-I 和 PA-II）。它们与感知机类似，不需要学习率。然而，与感知机不同的是，它们包含一个正则化参数 eta0（参考论文中的 \(C\)）。

对于分类，SGDClassifier(loss="hinge", penalty=None, learning_rate="pa1", eta0=1.0) 可用于 PA-I，或者使用 learning_rate="pa2" 用于 PA-II。对于回归，SGDRegressor(loss="epsilon_insensitive", penalty=None, learning_rate="pa1", eta0=1.0) 可用于 PA-I，或者使用 learning_rate="pa2" 用于 PA-II。

参考文献...

• “Online Passive-Aggressive Algorithms” K. Crammer, O. Dekel, J. Keshat, S. Shalev-Shwartz, Y. Singer - JMLR 7 (2006)

1.1.14. 稳健回归：异常值和建模错误

稳健回归旨在在存在损坏数据（异常值或模型错误）的情况下拟合回归模型。

1.1.14.1. 不同场景和有用概念

处理被异常值损坏的数据时，需要注意以下几点：

• X 或 y 中的异常值?

  y 方向上的异常值	X 方向上的异常值

• 异常值的比例与误差幅度

  异常点的数量很重要，异常的程度也很重要。

  小异常值	大异常值

稳健拟合的一个重要概念是崩溃点 (breakdown point)：指在拟合开始偏离正常数据之前，数据集中可以包含异常值的比例。

请注意，通常在高维设置（大 n_features）下进行稳健拟合非常困难。此处的稳健模型在这些设置下可能无法发挥作用。

权衡：选择哪个估计器？

Scikit-learn 提供了 3 种稳健回归估计器：RANSAC、Theil Sen 和 HuberRegressor。

• 除非样本数量非常大（即 n_samples >> n_features），否则 HuberRegressor 应该比 RANSAC 和 Theil Sen 更快。这是因为 RANSAC 和 Theil Sen 是在数据的较小子集上进行拟合的。然而，在使用默认参数时，Theil Sen 和 RANSAC 可能不如 HuberRegressor 稳健。

• RANSAC 比 Theil Sen 更快，并且在处理样本数量增加时扩展性更好。

• RANSAC 可以更好地处理 y 方向上的大异常值（这是最常见的情况）。

• Theil Sen 可以更好地处理 X 方向上的中等大小异常值，但这种特性在高维设置下会消失。

如有疑问，请使用 RANSAC。

1.1.14.2. RANSAC：随机采样一致性 (RANdom SAmple Consensus)

RANSAC (随机采样一致性) 从完整数据集的随机内点子集中拟合模型。

RANSAC 是一种非确定性算法，仅以一定概率产生合理结果，该概率取决于迭代次数（请参阅 max_trials 参数）。它通常用于线性和非线性回归问题，在摄影测量计算机视觉领域尤其受欢迎。

该算法将完整的输入样本数据分成一组内点（可能包含噪声）和异常值（例如由错误测量或对数据的无效假设引起）。最终的模型仅根据确定的内点进行估计。

示例

• 使用 RANSAC 进行鲁棒线性模型估计

• 鲁棒线性估计器拟合

算法细节...

每次迭代执行以下步骤：

1. 从原始数据中选择 min_samples 个随机样本，并检查该数据子集是否有效（请参阅 is_data_valid）。

2. 将模型拟合到随机子集（estimator.fit）并检查估计的模型是否有效（请参阅 is_model_valid）。

3. 通过计算估计模型的残差（estimator.predict(X) - y）将所有数据分类为内点或异常值——所有绝对残差小于或等于 residual_threshold 的数据样本均被视为内点。

4. 如果内点样本数量最大，则将拟合的模型保存为最佳模型。如果当前估计的模型具有相同数量的内点，则仅当其得分更好时才将其视为最佳模型。

这些步骤执行最大次数（max_trials）或直到满足特殊停止条件之一（请参阅 stop_n_inliers 和 stop_score）。最终模型使用之前确定的最佳模型的所有内点样本（一致集）进行估计。

is_data_valid 和 is_model_valid 函数允许识别并拒绝退化的随机子样本组合。如果不需要估计的模型来识别退化情况，则应使用 is_data_valid，因为它在拟合模型之前被调用，从而带来更好的计算性能。

参考文献...

• https://en.wikipedia.org/wiki/RANSAC

• “Random Sample Consensus: A Paradigm for Model Fitting with Applications to Image Analysis and Automated Cartography” Martin A. Fischler and Robert C. Bolles - SRI International (1981)

• “Performance Evaluation of RANSAC Family” Sunglok Choi, Taemin Kim and Wonpil Yu - BMVC (2009)

1.1.14.3. Theil-Sen 估计器：基于广义中位数的估计器

TheilSenRegressor 估计器使用了多维中位数的推广，因此对多变量异常值具有稳健性。但请注意，估计器的稳健性会随着问题维数的增加而迅速降低。在高维情况下，它会丧失稳健性，性能不优于普通最小二乘法。

示例

• Theil-Sen 回归

• 鲁棒线性估计器拟合

理论考量...

TheilSenRegressor 在渐近效率和无偏估计方面与 普通最小二乘法 (OLS) 相当。与 OLS 不同，Theil-Sen 是一种非参数方法，这意味着它不对数据的底层分布做任何假设。由于 Theil-Sen 是一种基于中位数的估计器，它对损坏数据（即异常值）更加稳健。在单变量设置中，Theil-Sen 在简单线性回归情况下的崩溃点约为 29.3%，这意味着它可以容忍高达 29.3% 的任意损坏数据。

scikit-learn 中 TheilSenRegressor 的实现遵循了多变量线性回归模型的推广 [14]，该模型使用了空间中位数，即中位数在多维空间中的推广 [15]。

在时间和空间复杂度方面，Theil-Sen 的扩展方式为：

\[\binom{n_{\text{samples}}}{n_{\text{subsamples}}}\]

这使得它在样本和特征数量众多的问题上无法穷举应用。因此，可以选择一个子群体大小，通过仅考虑所有可能组合的一个随机子集来限制时间和空间复杂度。

References

[14]

Xin Dang, Hanxiang Peng, Xueqin Wang and Heping Zhang: Theil-Sen Estimators in a Multiple Linear Regression Model.

[15]

20. Kärkkäinen and S. Äyrämö: On Computation of Spatial Median for Robust Data Mining.

另请参阅 Wikipedia 页面

1.1.14.4. Huber 回归

HuberRegressor 与 Ridge 不同，因为它对由 epsilon 参数定义为异常值的样本应用线性损失。如果样本的绝对误差小于阈值 epsilon，则该样本被归类为内点。它与 TheilSenRegressor 和 RANSACRegressor 不同，因为它不会忽略异常值的影响，而是赋予其较小的权重。

示例

• HuberRegressor vs Ridge 在具有强离群点的数据集上的比较

数学细节...

HuberRegressor 最小化：

\[\min_{w, \sigma} {\sum_{i=1}^n\left(\sigma + H_{\epsilon}\left(\frac{X_{i}w - y_{i}}{\sigma}\right)\sigma\right) + \alpha {||w||_2}^2}\]

其中损失函数为：

\[\begin{split}H_{\epsilon}(z) = \begin{cases} z^2, & \text {if } |z| < \epsilon, \\ 2\epsilon|z| - \epsilon^2, & \text{otherwise} \end{cases}\end{split}\]

建议将 epsilon 参数设置为 1.35 以达到 95% 的统计效率。

References

• Peter J. Huber, Elvezio M. Ronchetti: Robust Statistics, Concomitant scale estimates, p. 172.

HuberRegressor 与将 SGDRegressor 的损失函数设置为 huber 的方式在以下方面有所不同：

• HuberRegressor 具有缩放不变性。一旦设置了 epsilon，将 X 和 y 缩放不同的值将产生与之前相同的异常值稳健性。相比之下，SGDRegressor 在 X 和 y 被缩放时，必须重新设置 epsilon。

• HuberRegressor 在样本数量较少的数据上使用效率更高，而 SGDRegressor 需要多次遍历训练数据才能达到相同的稳健性。

请注意，此估计器与 R 语言的稳健回归实现 不同，因为 R 实现执行的是加权最小二乘法，其权重是根据每个样本的残差是否超过特定阈值来赋予的。

1.1.15. 分位数回归 (Quantile Regression)

分位数回归在给定 \(X\) 的条件下估计 \(y\) 的中位数或其他分位数，而普通最小二乘法 (OLS) 估计的是条件均值。

如果您对预测区间而不是点预测感兴趣，分位数回归可能很有用。有时，预测区间是基于预测误差服从均值为零且方差恒定的正态分布的假设计算的。即使对于具有非恒定（但可预测）方差或非正态分布的误差，分位数回归也能提供合理的预测区间。

基于最小化 pinball 损失，条件分位数也可以由线性模型以外的模型进行估计。例如，如果将 GradientBoostingRegressor 的 loss 参数设置为 "quantile"，并将 alpha 参数设置为要预测的分位数，则它可以预测条件分位数。请参阅 梯度提升回归的预测区间 中的示例。

大多数分位数回归实现都基于线性规划问题。当前的实现基于 scipy.optimize.linprog。

示例

• 分位数回归

数学细节...

作为一个线性模型，QuantileRegressor 给出了第 \(q\) 个分位数的线性预测 \(\hat{y}(w, X) = Xw\)，其中 \(q \in (0, 1)\)。权重或系数 \(w\) 然后通过以下最小化问题找到：

\[\min_{w} {\frac{1}{n_{\text{samples}}} \sum_i PB_q(y_i - X_i w) + \alpha ||w||_1}.\]

这包括 pinball 损失（也称为线性损失），另请参阅 mean_pinball_loss：

\[\begin{split}PB_q(t) = q \max(t, 0) + (1 - q) \max(-t, 0) = \begin{cases} q t, & t > 0, \\ 0, & t = 0, \\ (q-1) t, & t < 0 \end{cases}\end{split}\]

以及由 alpha 参数控制的 L1 惩罚项，类似于 Lasso。

由于 pinball 损失在残差中只是线性的，因此分位数回归比基于平方误差的均值估计对异常值更加稳健。处于两者之间的是 HuberRegressor。

参考文献...

• Koenker, R., & Bassett Jr, G. (1978). Regression quantiles. Econometrica: journal of the Econometric Society, 33-50.

• Portnoy, S., & Koenker, R. (1997). The Gaussian hare and the Laplacian tortoise: computability of squared-error versus absolute-error estimators. Statistical Science, 12, 279-300.

• Koenker, R. (2005). Quantile Regression. Cambridge University Press.

1.1.16. 多项式回归：利用基函数扩展线性模型

机器学习中的一种常见模式是使用在数据的非线性函数上训练的线性模型。这种方法保持了线性方法通常较快的性能，同时允许它们拟合范围更广的数据。

数学细节...

例如，简单的线性回归可以通过从系数构造多项式特征来扩展。在标准线性回归情况下，对于二维数据，模型可能如下所示：

\[\hat{y}(w, x) = w_0 + w_1 x_1 + w_2 x_2\]

如果我们想拟合一个抛物面而不是平面，我们可以将特征组合成二阶多项式，使模型看起来如下所示：

\[\hat{y}(w, x) = w_0 + w_1 x_1 + w_2 x_2 + w_3 x_1 x_2 + w_4 x_1^2 + w_5 x_2^2\]

（有时令人惊讶的）观察结果是，这仍然是一个线性模型：为了看到这一点，想象一下创建一组新的特征：

\[z = [x_1, x_2, x_1 x_2, x_1^2, x_2^2]\]

通过这种数据的重新标记，我们的问题可以写成：

\[\hat{y}(w, z) = w_0 + w_1 z_1 + w_2 z_2 + w_3 z_3 + w_4 z_4 + w_5 z_5\]

我们可以看到，结果的多项式回归与我们上面考虑的线性模型属于同一类别（即模型在 \(w\) 中是线性的），并且可以用相同的技术解决。通过在由这些基函数构建的更高维空间中考虑线性拟合，模型具有拟合更广泛数据的灵活性。

这是将此想法应用于一维数据，使用不同次数的多项式特征的示例：

此图是使用 PolynomialFeatures 转换器创建的，它将输入数据矩阵转换为给定次数的新数据矩阵。它可以按如下方式使用：

>>> from sklearn.preprocessing import PolynomialFeatures
>>> import numpy as np
>>> X = np.arange(6).reshape(3, 2)
>>> X
array([[0, 1],
       [2, 3],
       [4, 5]])
>>> poly = PolynomialFeatures(degree=2)
>>> poly.fit_transform(X)
array([[ 1.,  0.,  1.,  0.,  0.,  1.],
       [ 1.,  2.,  3.,  4.,  6.,  9.],
       [ 1.,  4.,  5., 16., 20., 25.]])

X 的特征已从 \([x_1, x_2]\) 转换为 \([1, x_1, x_2, x_1^2, x_1 x_2, x_2^2]\)，现在可以在任何线性模型中使用。

这种预处理可以使用 Pipeline 工具进行简化。可以创建一个表示简单多项式回归的单个对象，并按如下方式使用：

>>> from sklearn.preprocessing import PolynomialFeatures
>>> from sklearn.linear_model import LinearRegression
>>> from sklearn.pipeline import Pipeline
>>> import numpy as np
>>> model = Pipeline([('poly', PolynomialFeatures(degree=3)),
...                   ('linear', LinearRegression(fit_intercept=False))])
>>> # fit to an order-3 polynomial data
>>> x = np.arange(5)
>>> y = 3 - 2 * x + x ** 2 - x ** 3
>>> model = model.fit(x[:, np.newaxis], y)
>>> model.named_steps['linear'].coef_
array([ 3., -2.,  1., -1.])

在多项式特征上训练的线性模型能够精确恢复输入的多项式系数。

在某些情况下，不需要包含任何单个特征的更高幂，而只需要将至多 \(d\) 个不同特征相乘的所谓交互特征。这些可以从 PolynomialFeatures 设置 interaction_only=True 中获得。

例如，在处理布尔特征时，对于所有 \(n\)，\(x_i^n = x_i\)，因此是无用的；但 \(x_i x_j\) 表示两个布尔值的合取。这样，我们可以用线性分类器解决 XOR 问题：

>>> from sklearn.linear_model import Perceptron
>>> from sklearn.preprocessing import PolynomialFeatures
>>> import numpy as np
>>> X = np.array([[0, 0], [0, 1], [1, 0], [1, 1]])
>>> y = X[:, 0] ^ X[:, 1]
>>> y
array([0, 1, 1, 0])
>>> X = PolynomialFeatures(interaction_only=True).fit_transform(X).astype(int)
>>> X
array([[1, 0, 0, 0],
       [1, 0, 1, 0],
       [1, 1, 0, 0],
       [1, 1, 1, 1]])
>>> clf = Perceptron(fit_intercept=False, max_iter=10, tol=None,
...                  shuffle=False).fit(X, y)

分类器的“预测”是完美的。

>>> clf.predict(X)
array([0, 1, 1, 0])
>>> clf.score(X, y)
1.0