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\]
../_images/sphx_glr_plot_ols_ridge_001.png

LinearRegression 在其 fit 方法中接受参数 Xysample_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. Ridge 回归和分类#

1.1.2.1. 回归#

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

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

复杂度参数 \(\alpha \geq 0\) 控制收缩量:\(\alpha\) 值越大,收缩量越大,从而系数对共线性越鲁棒。

../_images/sphx_glr_plot_ridge_path_001.png

与其它线性模型一样,Ridge 在其 fit 方法中接受数组 Xy,并将其线性模型的系数 \(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’

以上条件均不满足。

示例

1.1.2.2. 分类#

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

使用(惩罚)最小二乘损失来拟合分类模型而不是更传统的逻辑损失或铰链损失似乎有些可疑。然而,在实践中,所有这些模型在准确性或精确率/召回率方面都可以得到相似的交叉验证分数,而 RidgeClassifier 使用的惩罚最小二乘损失允许选择不同的数值求解器,具有不同的计算性能特征。

RidgeClassifier 的速度可能比例如 LogisticRegression 快得多,尤其是当类别数量很多时,因为它只需计算一次投影矩阵 \((X^T X)^{-1} X^T\)

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

示例

1.1.2.3. Ridge 复杂度#

此方法的时间复杂度与 普通最小二乘法 的复杂度相同。

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

RidgeCVRidgeClassifierCV 实现了带内置 \(alpha\) 参数交叉验证的 Ridge 回归/分类。它们的工作方式与 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]], [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.01)

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

参考文献#
[RL2007]

“Notes on Regularized Least Squares”,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\) 路径上的系数。

示例

注意

Lasso 特征选择

由于 Lasso 回归会产生稀疏模型,因此它可以用于特征选择,如 基于 L1 的特征选择 中所述。

参考文献#

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

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\)。间隔安全筛选规则就是这样的一个工具。在优化算法的任何时候,它们都可以告诉我们哪些特征可以安全地排除,即确定地设置为零。

参考文献#

1.1.3.2. 设置正则化参数#

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

1.1.3.2.1. 使用交叉验证#

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

对于具有许多共线性特征的高维数据集,LassoCV 通常更受青睐。然而,LassoLarsCV 在探索 alpha 参数的更多相关值方面具有优势,并且如果样本数量远小于特征数量,它通常比 LassoCV 更快。

lasso_cv_1 lasso_cv_2

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

或者,估计器 LassoLarsIC 提出了使用赤池信息准则(AIC)和贝叶斯信息准则(BIC)。它是一种计算成本更低的查找最佳 \(alpha\) 值的方法,因为正则化路径仅计算一次,而不是使用 k 折交叉验证时的 k+1 次。

实际上,这些准则是在样本内训练集上计算的。简而言之,它们通过模型的灵活性来惩罚不同 Lasso 模型过于乐观的分数(参见下面的“数学细节”部分)。

然而,这类准则需要正确估计解的自由度,是针对大样本(渐近结果)推导出来的,并假设正在调查的模型是候选模型。当问题病态时(例如,特征数量多于样本数量),它们也往往会失效。

../_images/sphx_glr_plot_lasso_lars_ic_001.png

示例

1.1.3.2.3. AIC 和 BIC 准则#

AIC(因此 BIC)的定义在文献中可能有所不同。在本节中,我们将提供有关 scikit-learn 中计算的准则的更多信息。

数学细节#

AIC 准则定义为

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

其中 \(\hat{L}\) 是模型的最大似然度,\(d\) 是参数的数量(也称为前述的自由度)。

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

\[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

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

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

1.1.4. 多任务 Lasso#

MultiTaskLasso 是一个线性模型,它联合估计多个回归问题的稀疏系数:y 是一个二维数组,形状为 (n_samples, n_tasks)。约束条件是,选定的特征对于所有回归问题(也称为任务)都是相同的。

下图比较了使用简单 Lasso 或 MultiTaskLasso 得到的系数矩阵 W 中非零项的位置。Lasso 估计产生分散的非零项,而 MultiTaskLasso 的非零项是完整的列。

multi_task_lasso_1 multi_task_lasso_2

拟合时间序列模型,强制要求任何激活的特征在所有时间都激活。

示例

数学细节#

在数学上,它由一个线性模型组成,该模型使用混合 \(\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. 弹性网络#

ElasticNet 是一个线性回归模型,它同时使用 \(\ell_1\)\(\ell_2\) 范数正则化系数进行训练。这种组合允许学习一个稀疏模型,其中只有少数权重是非零的,就像 Lasso 一样,同时仍然保持 Ridge 的正则化特性。我们通过 l1_ratio 参数控制 \(\ell_1\)\(\ell_2\) 的凸组合。

当存在多个相互关联的特征时,弹性网络很有用。Lasso 倾向于随机选择其中一个,而弹性网络则倾向于同时选择两者。

在 Lasso 和 Ridge 之间进行权衡的一个实际优势是,它允许弹性网络继承 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}\]
../_images/sphx_glr_plot_lasso_lasso_lars_elasticnet_path_002.png

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

示例

参考文献#

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

  • “Regularization Path For Generalized linear Models by Coordinate Descent”,Friedman, Hastie & Tibshirani, J Stat Softw, 2010(论文)。

  • “An Interior-Point Method for Large-Scale L1-Regularized Least Squares,” S. J. Kim, K. Koh, M. Lustig, S. Boyd and D. Gorinevsky, in IEEE Journal of Selected Topics in Signal Processing, 2007(论文

1.1.6. 多任务弹性网络#

MultiTaskElasticNet 是一个弹性网络模型,它联合估计多个回归问题的稀疏系数:Y 是一个形状为 (n_samples, n_tasks) 的二维数组。约束条件是,选定的特征对于所有回归问题(也称为任务)都是相同的。

在数学上,它由一个线性模型组成,该模型使用混合 \(\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_pathlars_path_gram 来使用。

1.1.8. LARS Lasso#

LassoLars 是一个使用 LARS 算法实现的 Lasso 模型,与基于坐标下降的实现不同,它产生精确解,该解是其系数范数的函数的分段线性曲线。

../_images/sphx_glr_plot_lasso_lasso_lars_elasticnet_path_001.png
>>> 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.        ])

示例

LARS 算法几乎免费提供正则化参数的完整系数路径,因此常见的操作是使用函数 lars_pathlars_path_gram 来检索路径。

数学公式#

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

LARS 解不是返回一个向量,而是由一条曲线表示,该曲线表示参数向量的 \(\ell_1\) 范数对应的解。完整的系数路径存储在形状为 (n_features, max_features + 1)coef_path_ 数组中。第一列始终为零。

References

1.1.9. 正交匹配追踪(OMP)#

OrthogonalMatchingPursuitorthogonal_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)方法,但更好,因为在每次迭代中,使用对先前选择的字典元素空间的正交投影来重新计算残差。

示例

参考文献#

1.1.10. 贝叶斯回归#

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

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

为了获得完全概率化的模型,输出 \(y\) 被假定为围绕 \(X w\) 服从高斯分布

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

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

贝叶斯回归的优点是:

  • 它适应当前数据。

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

贝叶斯回归的缺点包括:

  • 模型推断可能耗时。

参考文献#

1.1.10.1. 贝叶斯 Ridge 回归#

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

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

\(\alpha\)\(\lambda\) 的先验被选择为伽马分布,这是高斯精度分布的共轭先验。所得模型称为贝叶斯 Ridge 回归,与经典 Ridge 类似。

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

此外还有四个超参数,即 \(\alpha\)\(\lambda\) 的伽马先验分布的 \(\alpha_1\)\(\alpha_2\)\(\lambda_1\)\(\lambda_2\)。这些通常选择为无信息的。默认值为 \(\alpha_1 = \alpha_2 = \lambda_1 = \lambda_2 = 10^{-6}\)

贝叶斯 Ridge 回归用于回归。

>>> 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])

由于贝叶斯框架,找到的权重与 普通最小二乘法 找到的权重略有不同。然而,贝叶斯 Ridge 回归对病态问题更鲁棒。

示例

参考文献#

1.1.10.2. 自动相关性确定 - ARD#

自动相关性确定(如 ARDRegression 中实现的)是一种线性模型,它与 贝叶斯 Ridge 回归 非常相似,但它能产生更稀疏的系数 \(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}\}\)

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

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

请参阅 比较线性贝叶斯回归模型 以了解 ARD 与 贝叶斯 Ridge 回归 的详细比较。

请参阅 稀疏信号的 L1 模型 以了解 Lasso、ARD 和 ElasticNet 等方法在相关数据上的比较。

References

1.1.11. 逻辑回归#

逻辑回归在 LogisticRegression 中实现。尽管有其名称,但根据 scikit-learn/ML 的命名法,它被实现为用于分类的线性模型,而不是回归。逻辑回归在文献中也称为 logit 回归、最大熵分类(MaxEnt)或对数线性分类器。在此模型中,单个试验的可能结果的概率使用逻辑函数进行建模。

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

注意

正则化

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

注意

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

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

示例

1.1.11.1. 二元情况#

为方便表示,我们假设对于数据点 \(i\),目标 \(y_i\) 的值为 \(\{0, 1\}\)。拟合后,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\)

我们目前通过 Cl1_ratio 参数提供了四种正则化或惩罚项 \(r(w)\) 的选择:

惩罚

\(r(w)\)

无(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\)

弹性网络(0<l1_ratio<1

\(\frac{1 - \rho}{2}w^T w + \rho \|w\|_1\)

对于弹性网络,\(\rho\)(对应于 l1_ratio 参数)控制 \(\ell_1\) 正则化与 \(\ell_2\) 正则化的强度。\(\rho = 1\) 时弹性网络等价于 \(\ell_1\)\(\rho=0\) 时等价于 \(\ell_2\)

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

1.1.11.2. 多项情况#

二元情况可以扩展到 \(K\) 个类别,得到多项逻辑回归,也可参阅 对数线性模型

注意

利用所有类别概率之和必须为一的事实,可以仅使用 \((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}\)

我们目前通过 Cl1_ratio 参数提供了四种正则化或惩罚项 \(r(W)\) 的选择,其中 \(m\) 是特征的数量。

惩罚

\(r(W)\)

无(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\)

弹性网络(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 惩罚

弹性网络(L1 + L2)

无惩罚

多类别支持

多项式多类

行为

惩罚截距(不良)

对于大型数据集更快

对未缩放的数据集鲁棒

“lbfgs”求解器默认使用,因为它具有鲁棒性。对于 n_samples >> n_features,“newton-cholesky”是一个不错的选择,并且可以达到很高的精度(tol 值很小)。对于大型数据集,“saga”求解器通常比“lbfgs”更快(尤其是在低精度(高 tol)时)。对于大型数据集,您也可以考虑使用 SGDClassifier 配合 loss="log_loss",这可能更快,但需要更多调整。

1.1.11.3.1. 求解器之间的差异#

fit_intercept=False 并且拟合的 coef_(或)要预测的数据为零时,LogisticRegressionLinearSVC 与外部 liblinear 库直接使用 solver=liblinear 时,得分可能会有所不同。这是因为对于 decision_function 为零的样本,LogisticRegressionLinearSVC 预测负类别,而 liblinear 预测正类别。请注意,具有 fit_intercept=False 且有许多 decision_function 为零的样本的模型,很可能是一个欠拟合、糟糕的模型,建议您将 fit_intercept=True 并增加 intercept_scaling

求解器详情#
  • 求解器 “liblinear” 使用坐标下降(CD)算法,并依赖于优秀的 C++ LIBLINEAR 库,该库随 scikit-learn 一起提供。但是,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)。因此,它是稀疏多项逻辑回归的首选求解器。它也是唯一支持弹性网络(0 < l1_ratio < 1)的求解器。

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

  • “newton-cholesky”求解器是一个精确的牛顿求解器,它计算 Hessian 矩阵并求解由此产生的线性系统。它对于 n_samples >> n_features 来说是一个很好的选择,并且可以达到很高的精度(tol 值很小),但有几个缺点:仅支持 \(\ell_2\) 正则化。此外,由于 Hessian 矩阵是显式计算的,因此内存使用量与 n_featuresn_classes 成二次方关系。

有关这些求解器的一些比较,请参阅 [9]

References

注意

稀疏逻辑回归特征选择

带有 \(\ell_1\) 惩罚的逻辑回归会产生稀疏模型,因此可以用于特征选择,如 基于 L1 的特征选择 中所述。

注意

P 值估计

在没有惩罚的回归情况下,可以获取系数的 P 值和置信区间。statsmodels 包 原生支持此功能。在 sklearn 中,也可以使用自助法。

LogisticRegressionCV 实现了带内置交叉验证支持的逻辑回归,以根据 scoring 属性查找最佳 Cl1_ratio 参数。对于高维密集数据,“newton-cg”、“sag”、“saga”和“lbfgs”求解器通常更快,因为它们支持预热启动(请参阅 术语表)。

1.1.12. 广义线性模型#

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

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

其次,平方损失函数被指数族分布的单位偏差 \(d\)(或更精确地说,指数离散模型(EDM)[11] 的再生指数离散模型)替换。

最小化问题变为

\[\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)在下图所示,

../_images/poisson_gamma_tweedie_distributions.png

服从泊松、Tweedie(power=1.5)和伽马分布的随机变量 Y 的 PDF,具有不同的均值(\(\mu\))。请注意泊松分布和 Tweedie(power=1.5)分布在 \(Y=0\) 处存在点质量,而伽马分布则不存在,因为它具有严格正的目标域。#

伯努利分布是一个离散概率分布,用于模拟伯努利试验——即只有两个互斥结果的事件。多项式分布是伯努利分布对于分类随机变量的推广。虽然伯努利分布中的随机变量有两个可能的结果,但多项式随机变量可以取 K 个可能类别中的一个,每个类别的概率是单独指定的。

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

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

  • 如果目标值为正值且偏斜,可以尝试使用带有对数链接的伽马分布。

  • 如果目标值比伽马分布的尾部更重,可以尝试使用逆高斯分布(甚至 Tweedie 族更高方差的幂)。

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

用例#
  • 农业/天气建模:每年降雨事件的数量(泊松),每次降雨事件的量(伽马),每年总降雨量(Tweedie / 复合泊松伽马)。

  • 风险建模/保险定价:每年索赔事件/保单持有人数量(泊松),每次事件成本(伽马),每位保单持有人每年总成本(Tweedie / 复合泊松伽马)。

  • 信用违约:贷款无法偿还的概率(伯努利)。

  • 欺诈检测:诸如现金转账等金融交易是欺诈交易的概率(伯努利)。

  • 预测性维护:每年生产中断事件的数量(泊松),中断持续时间(伽马),每年总中断时间(Tweedie / 复合泊松伽马)。

  • 药物测试:在一组试验中治愈患者的概率,或患者出现副作用的概率(伯努利)。

  • 新闻分类:将新闻文章分类到三个类别,即财经新闻、政治新闻和娱乐新闻(多项式)。

References

1.1.12.1. 用法#

TweedieRegressor 实现了一个用于 Tweedie 分布的广义线性模型,该模型允许使用适当的 power 参数来建模上述任何分布。特别是

  • power = 0:正态分布。在此情况下,RidgeElasticNet 等特定估计器通常更合适。

  • 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)

示例

实际注意事项#

特征矩阵 X 在拟合之前应进行标准化。这可确保惩罚项平等对待特征。

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

如果要建模相对频率,即每暴露量(时间、体积……)的计数,可以通过使用泊松分布并将 \(y=\frac{\mathrm{counts}}{\mathrm{exposure}}\) 作为目标值,同时将 \(\mathrm{exposure}\) 作为样本权重来完成。有关具体示例,请参见例如 Tweedie 回归在保险索赔中的应用

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

1.1.13. 随机梯度下降 - SGD#

随机梯度下降是一种简单但非常有效的拟合线性模型的方法。当样本数量(以及特征数量)非常大时,它特别有用。partial_fit 方法允许在线/离核学习。

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

有关更多详细信息,您可以参考专用的 随机梯度下降 文档部分。

1.1.13.1. 感知器#

Perceptron 是另一种简单的分类算法,适用于大规模学习,并源自 SGD。默认情况下

  • 它不需要学习率。

  • 它没有进行正则化(惩罚)。

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

最后一个特性意味着感知器比具有 hinge 损失的 SGD 训练速度稍快,并且生成的模型更稀疏。

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

1.1.13.2. 被动攻击算法#

被动攻击(PA)算法是另一种(PA-I 和 PA-II)用于大规模在线学习的算法族,它源自 SGD。它们与感知器相似,因为它们不需要学习率。然而,与感知器不同的是,它们包含一个正则化参数 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. 鲁棒回归:异常值和建模错误#

鲁棒回归旨在在数据被污染(存在异常值或模型错误)的情况下拟合回归模型。

../_images/sphx_glr_plot_theilsen_001.png

1.1.14.1. 不同场景和实用概念#

在处理被异常值污染的数据时,需要考虑几点:

  • X 中的异常值或 y 中的异常值?

    y 方向的异常值

    X 方向的异常值

    y_outliers

    X_outliers

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

    异常点的数量很重要,但它们异常的程度也很重要。

    小异常值

    大异常值

    y_outliers

    large_y_outliers

鲁棒拟合的一个重要概念是“失效点”:允许多少比例的数据是异常的,以至于拟合开始偏离内点数据。

请注意,通常情况下,在高维设置(n_features 很大)下进行鲁棒拟合非常困难。在此类设置中,这里的鲁棒模型可能无法正常工作。

1.1.14.2. RANSAC:随机样本共识#

RANSAC(RANdom SAmple Consensus)从完整数据集的随机内点子集中拟合模型。

RANSAC 是一种非确定性算法,只能以一定的概率产生合理的结果,该概率取决于迭代次数(参见 max_trials 参数)。它通常用于线性和非线性回归问题,在摄影测量计算机视觉领域尤为流行。

该算法将完整的输入样本数据分为一组内点,这些内点可能存在噪声,以及一组异常点,例如由错误测量或无效数据假设引起。然后,仅根据确定的内点来估计结果模型。

../_images/sphx_glr_plot_ransac_001.png

示例

算法详情#

每次迭代执行以下步骤:

  1. 从原始数据中选择 min_samples 个随机样本,并检查数据集是否有效(参见 is_data_valid)。

  2. 使用随机子集拟合模型(estimator.fit),并检查估计的模型是否有效(参见 is_model_valid)。

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

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

这些步骤执行最多 max_trials 次,或直到满足特殊停止条件之一(参见 stop_n_inliersstop_score)。最终模型使用先前确定的最佳模型的所有内点样本(共识集)进行估计。

is_data_validis_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 估计器:基于广义中位数的估计器#

The TheilSenRegressor 估计器使用了多维中位数的推广。因此,它对多元异常值具有鲁棒性。但请注意,该估计器的鲁棒性会随着问题维度的增加而迅速下降。在高维情况下,它会失去其鲁棒性,并且不比普通最小二乘法好。

示例

理论考量#

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

../_images/sphx_glr_plot_theilsen_001.png

scikit-learn 中 TheilSenRegressor 的实现遵循了对多变量线性回归模型 [14] 的推广,使用了空间中位数,这是中位数向多维的推广 [15]

在时间和空间复杂度方面,Theil-Sen 的计算复杂度为

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

这使得它在处理大量样本和特征的问题时,穷举应用是不可行的。因此,可以通过仅考虑所有可能组合的随机子集来限制时间和空间复杂度。

References

另请参见 Wikipedia 页面

1.1.14.4. Huber 回归#

HuberRegressorRidge 的不同之处在于,它对由 epsilon 参数定义的异常值应用线性损失。如果样本的绝对误差小于阈值 epsilon,则该样本被归类为内点。它与 TheilSenRegressorRANSACRegressor 不同,因为它不忽略异常值的影响,而是赋予它们较低的权重。

../_images/sphx_glr_plot_huber_vs_ridge_001.png

示例

数学细节#

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 与使用 loss 设置为 huberSGDRegressor 的不同之处在于:

  • HuberRegressor 是尺度不变的。一旦设置了 epsilon,将 Xy 按不同值进行缩放,会产生与之前相同的鲁棒性。而对于 SGDRegressor,当 Xy 被缩放时,epsilon 需要重新设置。

  • HuberRegressor 在处理少量样本数据时应比 SGDRegressor 更高效,而 SGDRegressor 需要多次遍历训练数据才能产生相同的鲁棒性。

请注意,此估计器不同于 R 中的鲁棒回归实现,因为 R 实现采用了加权最小二乘法,根据残差大于特定阈值的程度为每个样本分配权重。

1.1.15. 分位数回归#

分位数回归估计 \(y\) 条件于 \(X\) 的中位数或其他分位数,而普通最小二乘法(OLS)估计条件均值。

如果人们对预测区间而不是点预测感兴趣,分位数回归可能很有用。有时,预测区间是基于预测误差服从零均值和恒定方差的正态分布的假设来计算的。即使对于误差具有非恒定(但可预测)方差或非正态分布的情况,分位数回归也能提供合理的预测区间。

../_images/sphx_glr_plot_quantile_regression_002.png

基于最小化 pinball loss,条件分位数也可以通过线性模型以外的模型进行估计。例如,GradientBoostingRegressor 可以预测条件分位数,如果其参数 loss 设置为 "quantile",并且参数 alpha 设置为要预测的分位数。有关示例,请参见 梯度提升回归的预测区间

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

示例

数学细节#

作为一个线性模型,QuantileRegressor\(q\)-th 分位数 \(q \in (0, 1)\) 提供线性预测 \(\hat{y}(w, X) = Xw\)。然后通过以下最小化问题找到权重或系数 \(w\)

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

这包括 pinball loss(也称为线性损失),另请参阅 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 loss 仅在线性残差中,因此分位数回归比基于平方误差的均值估计对异常值更鲁棒。介于两者之间的是 HuberRegressor

参考文献#

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\) 上是线性的),并且可以使用相同的技术来求解。通过在用这些基函数构建的更高维空间中考虑线性拟合,模型具有适应更广泛数据的灵活性。

这是一个将此思想应用于一维数据的示例,使用了不同次数的多项式特征:

../_images/sphx_glr_plot_polynomial_interpolation_001.png

该图是使用 PolynomialFeatures transformer 创建的,该 transformer 将输入数据矩阵转换为给定次数的新数据矩阵。它可如下使用:

>>> 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 来获得。

例如,当处理布尔特征时,\(x_i^n = x_i\) 对于所有 \(n\) 都成立,因此无用;但 \(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