线性回归

在这节课上我们将学习线性回归。 我们将先理解它背后的数学基础原理再用python去实现它。 我们还将通过方法去讲解线性模型。

概述

$\hat{y} = XW$

where:

  • $\hat{y}$ = 预测值 | $\in \mathbb{R}^{NX1}$ ($N$ 是样本的个数)
  • $X$ = 输入 | $\in \mathbb{R}^{NXD}$ ($D$ 是特征的个数)
  • $W$ = 权重 | $\in \mathbb{R}^{DX1}$
  • 目标: 通过线性模型的输入 $X$ 去预测 $\hat{y}$。模型将会寻找一条最优的线使得我们的预测值和目标值最为接近。训练数据 $(X, y)$ 用来训练这个模型并且通过随机梯度下降(SGD)学习权重 $W$。
  • 优点:
    • 计算简单。
    • 解释性强。
    • 可用于连续(continuous)和无序的类别(categorical)特征。
  • 缺点:
    • 线性模型只能用于线性可分的数据(针对于分类任务).
    • 但是通常来讲不会用于分类任务,仅仅用于回归问题。
  • 其他: 当然你也可以使用线性回归去做二分类任务,如果预测出的连续数值高于一个阈值它就属于一个特定的分类。但是我们在未来的课程中将会介绍可用于做二分类任务更好的模型,所以我们本次课程只会集中在怎么用线性回归去做回归任务。

训练

步骤:

  1. 随机初始化模型的权重$W$。
  2. 将输入值 $X$ 传入模型并且得到预测值$\hat{y}$。
  3. 通过损失函数来计算预测值$\hat{y}$和真实值$\hat{y}$之间的差距,从而得到损失值$J$。普遍在线性回归中用到的损失函数是均方误差(MSE)。这个函数计算出预测值和真实值之间的差距的平方($\frac{1}{2}$ 没有数学意义,只是在求导的时候可以正好和平方抵消,方便计算)。
    • $MSE = J(\theta) = \frac{1}{2}\sum_{i}(\hat{y}_i - y_i)^2$
  4. 计算出对于模型权重的损失梯度$J(\theta)$
    • $J(\theta) = \frac{1}{2}\sum_{i}(\hat{y}_i - y_i)^2 = \frac{1}{2}\sum_{i}(X_iW - y_i)^2 $
    • $\frac{\partial{J}}{\partial{W}} = X(\hat{y} - y)$
  5. 我们使用学习率$\alpha$和一个优化方法(比如随机梯度下降),通过反向传播来更新权重$W$。 一个简单的比方就是梯度可以告诉你在哪个方向上增加数值,然后通过减法来使得损失值$J(\theta)$越来越小。
    • $W = W- \alpha\frac{\partial{J}}{\partial{W}}$
  6. 重复2 - 4步直到模型表现最好(也可以说直到损失收敛)。

数据

我们将自己创建一些假数据应用在线性回归上。

In [1]:
from argparse import Namespace
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
In [2]:
# 参数
args = Namespace(
    seed=1234,
    data_file="sample_data.csv",
    num_samples=100,
    train_size=0.75,
    test_size=0.25,
    num_epochs=100,
)

# 设置随机种子来保证实验结果的可重复性。
np.random.seed(args.seed)
In [3]:
# 生成数据
def generate_data(num_samples):
    X = np.array(range(num_samples))
    y = 3.65*X + 10
    return X, y
In [4]:
# 生成随机数据
X, y = generate_data(args.num_samples)
data = np.vstack([X, y]).T
df = pd.DataFrame(data, columns=['X', 'y'])
df.head()
Out[4]:
X y
0 0.0 10.00
1 1.0 13.65
2 2.0 17.30
3 3.0 20.95
4 4.0 24.60
In [5]:
# 画散点图
plt.title("Generated data")
plt.scatter(x=df["X"], y=df["y"])
plt.show()

Scikit-learn 实现方法

注意: LinearRegression类在Scikit-learn中使用的是正规方程法来做的拟合。然而,我们将会使用Scikit-learn中的随机梯度下降SGDRegressor类来拟合数据。我们使用这个优化方法是因为在未来的几节课程中我们也会使用到它。

In [6]:
# 调包
from sklearn.linear_model.stochastic_gradient import SGDRegressor
from sklearn.preprocessing import StandardScaler
from sklearn.model_selection import train_test_split
In [7]:
# 划分数据到训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(
    df["X"].values.reshape(-1, 1), df["y"], test_size=args.test_size, 
    random_state=args.seed)
print ("X_train:", X_train.shape)
print ("y_train:", y_train.shape)
print ("X_test:", X_test.shape)
print ("y_test:", y_test.shape)
X_train: (75, 1)
y_train: (75,)
X_test: (25, 1)
y_test: (25,)

我们需要标准化我们的数据(零均值和单位方差),以便正确使用SGD并在速度上优化。

In [24]:
# 标准化训练集数据 (mean=0, std=1) 
X_scaler = StandardScaler().fit(X_train)
y_scaler = StandardScaler().fit(y_train.values.reshape(-1,1))

# 在训练集和测试集上进行标准化操作
standardized_X_train = X_scaler.transform(X_train)
standardized_y_train = y_scaler.transform(y_train.values.reshape(-1,1)).ravel()
standardized_X_test = X_scaler.transform(X_test)
standardized_y_test = y_scaler.transform(y_test.values.reshape(-1,1)).ravel()


# 检查
print ("mean:", np.mean(standardized_X_train, axis=0), 
       np.mean(standardized_y_train, axis=0)) # mean 应该是 ~0
print ("std:", np.std(standardized_X_train, axis=0), 
       np.std(standardized_y_train, axis=0))   # std 应该是 1
mean: [8.22952817e-17] -1.5617137213060536e-16
std: [1.] 0.9999999999999999
In [25]:
# 初始化模型
lm = SGDRegressor(loss="squared_loss", penalty="none", max_iter=args.num_epochs)
In [26]:
# 训练
lm.fit(X=standardized_X_train, y=standardized_y_train)
Out[26]:
SGDRegressor(alpha=0.0001, average=False, epsilon=0.1, eta0=0.01,
       fit_intercept=True, l1_ratio=0.15, learning_rate='invscaling',
       loss='squared_loss', max_iter=100, n_iter=None, penalty='none',
       power_t=0.25, random_state=None, shuffle=True, tol=None, verbose=0,
       warm_start=False)
In [27]:
# 预测 (还未标准化)
pred_train = (lm.predict(standardized_X_train) * np.sqrt(y_scaler.var_)) + y_scaler.mean_
pred_test = (lm.predict(standardized_X_test) * np.sqrt(y_scaler.var_)) + y_scaler.mean_

评估

有很多种方法可以来评估我们模型表现的好坏。

In [28]:
import matplotlib.pyplot as plt
In [29]:
# 训练和测试集上的均方误差 MSE
train_mse = np.mean((y_train - pred_train) ** 2)
test_mse = np.mean((y_test - pred_test) ** 2)
print ("train_MSE: {0:.2f}, test_MSE: {1:.2f}".format(train_mse, test_mse))
train_MSE: 0.00, test_MSE: 0.00

除了使用MSE,如果我们只有一个特征向量,我们可以把他们可视化出来直观的评估模型。

In [30]:
# 图例大小
plt.figure(figsize=(15,5))

# 画出训练数据
plt.subplot(1, 2, 1)
plt.title("Train")
plt.scatter(X_train, y_train, label="y_train")
plt.plot(X_train, pred_train, color="red", linewidth=1, linestyle="-", label="lm")
plt.legend(loc='lower right')

# 画出测试数据
plt.subplot(1, 2, 2)
plt.title("Test")
plt.scatter(X_test, y_test, label="y_test")
plt.plot(X_test, pred_test, color="red", linewidth=1, linestyle="-", label="lm")
plt.legend(loc='lower right')

# 显示图例
plt.show()

推论

In [31]:
# 传入我们自己的输入值
X_infer = np.array((0, 1, 2), dtype=np.float32)
standardized_X_infer = X_scaler.transform(X_infer.reshape(-1, 1))
pred_infer = (lm.predict(standardized_X_infer) * np.sqrt(y_scaler.var_)) + y_scaler.mean_
print (pred_infer)
df.head(3)
[10.00362276 13.65354177 17.30346078]
Out[31]:
X y
0 0.0 10.00
1 1.0 13.65
2 2.0 17.30

可解释性

线性回归有很强的可解释性。每一个特征都有一个系数来控制对输出值y的影响大小。我们可以这样解释这个系数: 如果我们把x增加1, 我们将把y增加 $W$ (~3.65)。

注意: 因为我们在梯度下降时标准化了我们的输入和输出,我们需要对我们的系数和截距做一个反标准化。过程可见下方。

In [32]:
# 未标准化系数
coef = lm.coef_ * (y_scaler.scale_/X_scaler.scale_)
intercept = lm.intercept_ * y_scaler.scale_ + y_scaler.mean_ - np.sum(coef*X_scaler.mean_)
print (coef) # ~3.65
print (intercept) # ~10
[3.64992205]
[10.00362489]

非标准化系数的证明:

注意我们的X和y都已经标准化了。

$\frac{\mathbb{E}[y] - \hat{y}}{\sigma_y} = W_0 + \sum_{j=1}^{k}W_jz_j$

$z_j = \frac{x_j - \bar{x}_j}{\sigma_j}$

$ \hat{y}_{scaled} = \frac{\hat{y}_{unscaled} - \bar{y}}{\sigma_y} = \hat{W_0} + \sum_{j=1}^{k} \hat{W}_j (\frac{x_j - \bar{x}_j}{\sigma_j}) $

$\hat{y}_{unscaled} = \hat{W}_0\sigma_y + \bar{y} - \sum_{j=1}^{k} \hat{W}_j(\frac{\sigma_y}{\sigma_j})\bar{x}_j + \sum_{j=1}^{k}(\frac{\sigma_y}{\sigma_j})x_j $

正则化

正规化有助于减少过拟合。下方是L2正则化(ridge regression)。有很多正则化的方法他们都可以使我们的模型减少过拟合。对于L2正则化, 我们会减小那些值很大的权重。数值很大的权重将会使模型更加看中它们的特征,但是我们希望的是模型会公平的对待所有的特征而不是仅仅权重很大的几个。 当然还有其他的正则化方法比如L1(lasso regression),它对于我们想创建更加稀疏的数据模型有好处,因为它会使得一些权重变成0,或者我们可以结合L2和L1正则化方法。

注意: 正则化不仅仅用于线性回归。它可以用于任何常规模型以及我们以后将会学到的模型。

  • $ J(\theta) = = \frac{1}{2}\sum_{i}(X_iW - y_i)^2 + \frac{\lambda}{2}\sum\sum W^2$
  • $ \frac{\partial{J}}{\partial{W}} = X (\hat{y} - y) + \lambda W $
  • $W = W- \alpha\frac{\partial{J}}{\partial{W}}$ where:
    • $\lambda$ 是正则化系数
In [33]:
# 初始化带有L2正则化的模型
lm = SGDRegressor(loss="squared_loss", penalty='l2', alpha=1e-2, 
                  max_iter=args.num_epochs)
In [34]:
# 训练
lm.fit(X=standardized_X_train, y=standardized_y_train)
Out[34]:
SGDRegressor(alpha=0.01, average=False, epsilon=0.1, eta0=0.01,
       fit_intercept=True, l1_ratio=0.15, learning_rate='invscaling',
       loss='squared_loss', max_iter=100, n_iter=None, penalty='l2',
       power_t=0.25, random_state=None, shuffle=True, tol=None, verbose=0,
       warm_start=False)
In [35]:
# 预测 (还未标准化)
pred_train = (lm.predict(standardized_X_train) * np.sqrt(y_scaler.var_)) + y_scaler.mean_
pred_test = (lm.predict(standardized_X_test) * np.sqrt(y_scaler.var_)) + y_scaler.mean_
In [36]:
# 训练集和测试集的MSE
train_mse = np.mean((y_train - pred_train) ** 2)
test_mse = np.mean((y_test - pred_test) ** 2)
print ("train_MSE: {0:.2f}, test_MSE: {1:.2f}".format(
    train_mse, test_mse))
train_MSE: 1.09, test_MSE: 1.15

正则化对于我们在做的这个数据帮助很少,因为我们在创建数据的时候用的就是一个线性的函数。但是对于现实中的数据,正则化就可以帮助我们构建更好的模型。

In [37]:
# 未标准化系数
coef = lm.coef_ * (y_scaler.scale_/X_scaler.scale_)
intercept = lm.intercept_ * y_scaler.scale_ + y_scaler.mean_ - (coef*X_scaler.mean_)
print (coef) # ~3.65
print (intercept) # ~10
[3.61384419]
[11.67785083]

类别变量

在我们的例子中,特征用的是连续的数值,那么假设我们要用类别的特征变量呢?一种选择就是使用独热编码来处理类别变量,这种方法用Pandas很容易实现,你可以用和上面相同的步骤来训练你的模型

In [38]:
# 创建类别特征
cat_data = pd.DataFrame(['a', 'b', 'c', 'a'], columns=['favorite_letter'])
cat_data.head()
Out[38]:
favorite_letter
0 a
1 b
2 c
3 a
In [39]:
dummy_cat_data = pd.get_dummies(cat_data) #独热编码 one-hot encoding,与dummy变量不同要注意。
dummy_cat_data.head()
Out[39]:
favorite_letter_a favorite_letter_b favorite_letter_c
0 1 0 0
1 0 1 0
2 0 0 1
3 1 0 0

现在你可以拼接上连续特征变量来训练线性模型。

TODO

  • 多项式回归
  • 一个简单的用正规方程的例子(sklearn.linear_model.LinearRegression)来分析优点和缺点,并且和随机梯度下降线性回归做对比。