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

一个普通最小二乘法简单和多重线性回归的实现。

> 从 mlxtend.regressor 导入线性回归

概述

简单线性回归模型的示例：

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

$$SSE = \sum_i \big(\text{目标}^{(i)} - \text{输出}^{(i)}\big)^2$$

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

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

正常方程（闭式解）

对于"较小"的数据集，应该优先考虑闭式解，因为在这种情况下计算（“昂贵的”）矩阵逆并不是一个问题。对于非常大的数据集，或者逆矩阵 $[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的梯度下降规则以获取详细信息。

随机洗牌的实现过程为：

参考文献

示例 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

# 可视化检查收敛性的成本并绘制线性模型：

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)

Ordinary least squares linear regression.

Parameters

Attributes

Examples

For usage examples, please see https://rasbt.github.io/mlxtend/user_guide/regressor/LinearRegression/

Methods

fit(X, y, init_params=True)

Learn model from training data.

Parameters

Returns

get_params(deep=True)

Get parameters for this estimator.

Parameters

Returns

predict(X)

Predict targets from X.

Parameters

Returns

set_params(params)

Set the parameters of this estimator. The method works on simple estimators as well as on nested objects (such as pipelines). The latter have parameters of the form <component>__<parameter> so that it's possible to update each component of a nested object.

Returns

self

adapted from https://github.com/scikit-learn/scikit-learn/blob/master/sklearn/base.py Author: Gael Varoquaux gael.varoquaux@normalesup.org License: BSD 3 clause

ython