LinearRegression: 普通最小二乘线性回归的实现

普通最小二乘简单线性回归和多元线性回归的实现。

from mlxtend.regressor import LinearRegression

概述

简单线性回归模型示例

在普通最小二乘 (OLS) 线性回归中，我们的目标是找到使垂直偏移最小的直线（或超平面）。换句话说，我们将最佳拟合线定义为使目标变量 (y) 和所有样本的预测输出之间的平方误差和 (SSE) 或均方误差 (MSE) 最小化的直线。 $i$ 在我们的数据集大小为 $n$ .

$SSE = \sum_i \big(\text{target}^{(i)} - \text{output}^{(i)}\big)^2$

$MSE = \frac{1}{n} \times SSE$

现在，LinearRegression 使用以下五种方法之一实现了用于执行普通最小二乘回归的线性回归模型：

正规方程
QR 分解法
SVD（奇异值分解）法
梯度下降
随机梯度下降

正规方程 (闭式解)

对于计算（昂贵的）矩阵逆不成问题的“较小”数据集，应首选闭式解。对于非常大的数据集，或逆矩阵 $[X^T X]$ 可能不存在（矩阵不可逆或奇异，例如在完全多重共线性情况下），应首选 QR、SVD 或梯度下降方法。

线性函数（线性回归模型）定义为

$y = w_0x_0 + w_1x_1 + ... + w_mx_m = \sum_{i=0}^{m} = \mathbf{w}^T\mathbf{x}$

其中 $y$ 是响应变量， $\mathbf{x}$ 是一个 $m$ 维样本向量，且 $\mathbf{w}$ 是权重向量（系数向量）。请注意， $w_0$ 表示模型的 y 轴截距，因此 $x_0=1$ .

使用闭式解（正规方程），我们按如下方式计算模型的权重

$\mathbf{w} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^Ty$

通过 QR 分解实现稳定的 OLS

QR 分解方法提供了一种比基于“正规方程”的闭式解析解数值更稳定的替代方案，并且可以更有效地计算大型矩阵的逆矩阵。

QR 分解法将给定矩阵分解为两个易于获得逆矩阵的矩阵。例如，给定矩阵 $X \in \mathbb{R}^{n \times m}$ ，其 QR 分解为两个矩阵是

$\mathbf{X = QR},$

其中

$Q \in \mathbb{R}^{n \times m}$ 是一个标准正交矩阵，使得 $Q^\top Q = QQ^\top = I$ 。第二个矩阵 $R \in \mathbf{R}^{m \times m}$ 是一个上三角矩阵。

然后可以按如下方式计算普通最小二乘回归模型的权重参数 [1]

$\mathbf{w} = \mathbf{R}^{-1}\mathbf{Q}^\top y.$

通过奇异值分解实现稳定的 OLS

另一种以数值稳定方式获取 OLS 模型权重的替代方法是奇异值分解 (SVD)，其定义为

$\mathbf{X}=\mathbf{U}\Sigma \mathbf{V}^\top,$

对于给定矩阵 $\mathbf{X}$ .

然后，可以证明 $X$ , $X^+$ 的伪逆可以按如下方式获得 [1]

$X^+ = \mathbf{U} \mathbf{\Sigma}^+ V^{\top}.$

请注意，虽然 $\Sigma$ 是包含 $\mathbf{X}$ , $\Sigma^{+}$ 的奇异值的对角矩阵，是包含奇异值倒数的对角矩阵。

然后可以按如下方式计算模型权重

$\mathbf{w} = \mathbf{X}^+ \mathbf{y}.$

请注意，这种 OLS 方法计算效率最低。然而，当直接方法（正规方程）或 QR 分解无法应用，或者正规方程（通过 $\mathbf{X}^T \mathbf{X}$ ）病态时 [3]，它是一种有用的方法。

梯度下降 (GD) 和随机梯度下降 (SGD)

有关详细信息，请参见梯度下降和随机梯度下降和推导线性回归和 Adaline 的梯度下降规则。

随机打乱的实现如下

对于一个或多个 epoch
- 随机打乱训练集中的样本
  - 对于训练样本 i
    - 计算梯度并执行权重更新

参考文献

[1] 第 3 章，第 55 页，回归的线性方法。Trevor Hastie; Robert Tibshirani; Jerome Friedman (2009). 统计学习基础：数据挖掘、推理和预测（第二版）。纽约：Springer。(ISBN 978–0–387–84858–7)
[2] G. Strang, 线性代数及其应用，第二版，佛罗里达州奥兰多，Academic Press, Inc.，1980 年，第 139-142 页。
[3] Douglas Wilhelm Harder. 工程数值分析。第 4.8 节，病态矩阵

示例 1 - 闭式解

import numpy as np
import matplotlib.pyplot as plt
from mlxtend.regressor import LinearRegression

X = np.array([ 1.0, 2.1, 3.6, 4.2, 6])[:, np.newaxis]
y = np.array([ 1.0, 2.0, 3.0, 4.0, 5.0])

ne_lr = LinearRegression()
ne_lr.fit(X, y)

print('Intercept: %.2f' % ne_lr.b_)
print('Slope: %.2f' % ne_lr.w_[0])

def lin_regplot(X, y, model):
    plt.scatter(X, y, c='blue')
    plt.plot(X, model.predict(X), color='red')    
    return

lin_regplot(X, y, ne_lr)
plt.show()

Intercept: 0.25
Slope: 0.81

png

示例 2 - QR 分解方法

import numpy as np
import matplotlib.pyplot as plt
from mlxtend.regressor import LinearRegression

X = np.array([ 1.0, 2.1, 3.6, 4.2, 6])[:, np.newaxis]
y = np.array([ 1.0, 2.0, 3.0, 4.0, 5.0])

qr_lr = LinearRegression(method='qr')
qr_lr.fit(X, y)

print('Intercept: %.2f' % qr_lr.b_)
print('Slope: %.2f' % qr_lr.w_[0])

def lin_regplot(X, y, model):
    plt.scatter(X, y, c='blue')
    plt.plot(X, model.predict(X), color='red')    
    return

lin_regplot(X, y, qr_lr)
plt.show()

Intercept: 0.25
Slope: 0.81

png

示例 3 - SVD 方法

import numpy as np
import matplotlib.pyplot as plt
from mlxtend.regressor import LinearRegression

X = np.array([ 1.0, 2.1, 3.6, 4.2, 6])[:, np.newaxis]
y = np.array([ 1.0, 2.0, 3.0, 4.0, 5.0])

svd_lr = LinearRegression(method='svd')
svd_lr.fit(X, y)

print('Intercept: %.2f' %svd_lr.b_)
print('Slope: %.2f' % svd_lr.w_[0])

def lin_regplot(X, y, model):
    plt.scatter(X, y, c='blue')
    plt.plot(X, model.predict(X), color='red')    
    return

lin_regplot(X, y, svd_lr)
plt.show()

Intercept: 0.25
Slope: 0.81

png

示例 4 - 梯度下降

import numpy as np
import matplotlib.pyplot as plt
from mlxtend.regressor import LinearRegression

X = np.array([ 1.0, 2.1, 3.6, 4.2, 6])[:, np.newaxis]
y = np.array([ 1.0, 2.0, 3.0, 4.0, 5.0])

gd_lr = LinearRegression(method='sgd',
                         eta=0.005, 
                         epochs=100,
                         minibatches=1,
                         random_seed=123,
                         print_progress=3)
gd_lr.fit(X, y)

print('Intercept: %.2f' % gd_lr.b_)
print('Slope: %.2f' % gd_lr.w_)

def lin_regplot(X, y, model):
    plt.scatter(X, y, c='blue')
    plt.plot(X, model.predict(X), color='red')    
    return

lin_regplot(X, y, gd_lr)
plt.show()

Iteration: 100/100 | Cost 0.08 | Elapsed: 0:00:00 | ETA: 0:00:00

Intercept: 0.22
Slope: 0.82

png

# Visualizing the cost to check for convergence and plotting the linear model:

plt.plot(range(1, gd_lr.epochs+1), gd_lr.cost_)
plt.xlabel('Epochs')
plt.ylabel('Cost')
plt.ylim([0, 0.2])
plt.tight_layout()
plt.show()

png

示例 5 - 随机梯度下降

import numpy as np
import matplotlib.pyplot as plt
from mlxtend.regressor import LinearRegression

X = np.array([ 1.0, 2.1, 3.6, 4.2, 6])[:, np.newaxis]
y = np.array([ 1.0, 2.0, 3.0, 4.0, 5.0])

sgd_lr = LinearRegression(method='sgd',
                          eta=0.01, 
                          epochs=100, 
                          random_seed=0, 
                          minibatches=len(y))
sgd_lr.fit(X, y)

print('Intercept: %.2f' % sgd_lr.w_)
print('Slope: %.2f' % sgd_lr.b_)

def lin_regplot(X, y, model):
    plt.scatter(X, y, c='blue')
    plt.plot(X, model.predict(X), color='red')    
    return

lin_regplot(X, y, sgd_lr)
plt.show()

Intercept: 0.82
Slope: 0.24

png

plt.plot(range(1, sgd_lr.epochs+1), sgd_lr.cost_)
plt.xlabel('Epochs')
plt.ylabel('Cost')
plt.ylim([0, 0.2])
plt.tight_layout()
plt.show()

png

示例 6 - 带小批量的随机梯度下降

import numpy as np
import matplotlib.pyplot as plt
from mlxtend.regressor import LinearRegression

X = np.array([ 1.0, 2.1, 3.6, 4.2, 6])[:, np.newaxis]
y = np.array([ 1.0, 2.0, 3.0, 4.0, 5.0])

sgd_lr = LinearRegression(method='sgd',
                          eta=0.01, 
                          epochs=100, 
                          random_seed=0, 
                          minibatches=3)
sgd_lr.fit(X, y)

print('Intercept: %.2f' % sgd_lr.b_)
print('Slope: %.2f' % sgd_lr.w_)

def lin_regplot(X, y, model):
    plt.scatter(X, y, c='blue')
    plt.plot(X, model.predict(X), color='red')    
    return

lin_regplot(X, y, sgd_lr)
plt.show()

Intercept: 0.24
Slope: 0.82

png

plt.plot(range(1, sgd_lr.epochs+1), sgd_lr.cost_)
plt.xlabel('Epochs')
plt.ylabel('Cost')
plt.ylim([0, 0.2])
plt.tight_layout()
plt.show()

png

API

LinearRegression(method='direct', eta=0.01, epochs=50, minibatches=None, random_seed=None, print_progress=0)

普通最小二乘线性回归。

参数

method : 字符串 (默认值: 'direct')

对于基于梯度下降的优化，使用 sgd（更多选项参见 minibatch 参数）。否则，如果为 direct（默认值），则使用解析法。对于替代的、数值更稳定的解决方案，可以使用 qr（QR 分解）或 svd（奇异值分解）。
eta : 浮点数 (默认值: 0.01)

求解器学习率（介于 0.0 和 1.0 之间）。与 method = 'sgd' 一起使用。（详细信息请参见 methods 参数）
epochs : 整数 (默认值: 50)

遍历训练数据集的次数。在每个 epoch 之前，如果 minibatches > 1，数据集会被打乱以防止随机梯度下降中的循环。与 method = 'sgd' 一起使用。（详细信息请参见 methods 参数）
minibatches : 整数 (默认值: None)

基于梯度优化的 mini-batch 数量。如果为 None：直接法、QR 或 SVD 法（详细信息参见 method 参数）如果为 1：梯度下降学习如果为 len(y)：随机梯度下降学习如果 1 < minibatches < len(y)：Mini-batch 学习
random_seed : 整数 (默认值: None)

设置用于打乱和初始化权重的随机状态。在 method = 'sgd' 中使用。（详细信息请参见 methods 参数）
print_progress : 整数 (默认值: 0)

如果 method = 'sgd'，则将拟合进度打印到 stderr。0：无输出 1：已用 epoch 数和成本 2：1 加已用时间 3：2 加预计完成时间

属性

w_ : 2维数组, shape={n_features, 1}

拟合后的模型权重。
b_ : 1维数组, shape={1,}

拟合后的偏置项。
cost_ : 列表

每个 epoch 后的平方误差和；如果 solver='normal equation' 则忽略

示例

使用示例请参见 https://mlxtend.cn/mlxtend/user_guide/regressor/LinearRegression/

方法

fit(X, y, init_params=True)

从训练数据中学习模型。

参数

X : {类数组, 稀疏矩阵}, shape = [n_samples, n_features]

训练向量，其中 n_samples 是样本数，n_features 是特征数。
y : 类数组, shape = [n_samples]

目标值。
init_params : 布尔值 (默认值: True)

在拟合之前重新初始化模型参数。设置为 False 以继续使用之前模型拟合的权重进行训练。

返回值

self : 对象

get_params(deep=True)

获取此估计器的参数。

参数

deep : 布尔值, 可选

如果为 True，则返回此估计器和包含的估计器子对象的参数。

返回值

params : 字符串到任意类型的映射

参数名称与其值之间的映射。

改编自 https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/base.py 作者：Gael Varoquaux gael.varoquaux@normalesup.org 许可：BSD 3 clause

predict(X)

从 X 预测目标。

参数

X : {类数组, 稀疏矩阵}, shape = [n_samples, n_features]

训练向量，其中 n_samples 是样本数，n_features 是特征数。

返回值

target_values : 类数组, shape = [n_samples]

预测的目标值。

set_params(params)

设置此估计器的参数。此方法适用于简单估计器以及嵌套对象（如管道）。后者具有 <组件>__<参数> 形式的参数，因此可以更新嵌套对象的每个组件。

返回值

self

改编自 https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/base.py 作者：Gael Varoquaux gael.varoquaux@normalesup.org 许可：BSD 3 clause

Python

按键	动作
`?`	打开帮助
`n`	下一页
`p`	上一页
`s`	搜索