44. 简单线性回归模型#

import numpy as np
import pandas as pd
import matplotlib as mpl
import matplotlib.pyplot as plt

FONTPATH = "fonts/SourceHanSerifSC-SemiBold.otf"
mpl.font_manager.fontManager.addfont(FONTPATH)
plt.rcParams['font.family'] = ['Source Han Serif SC']

以下的简单回归模型估计了两个变量 \(x_i\) 和 \(y_i\) 之间的关系

\[ y_i = \alpha + \beta x_i + \epsilon_i, i = 1,2,...,N \]

其中 \(\epsilon_i\) 表示最佳拟合线与给定 \(x_i\) 的样本值 \(y_i\) 之间的误差。

我们的目标是选择 \(\alpha\) 和 \(\beta\) 的值来为一些适用于变量 \(x_i\) 和 \(y_i\) 的数据构建一条“最佳”拟合线。

让我们假设存在一个简单的数据集，其中包含 10 个变量 \(x_i\) 和 \(y_i\) 的观测值：

\(y_i\)

\(x_i\)

1

2000

32

2

1000

21

3

1500

24

4

2500

35

5

500

10

6

900

11

7

1100

22

8

1500

21

9

1800

27

10

250

2

让我们把 \(y_i\) 视为一个冰淇淋车的销售额，而 \(x_i\) 是记录当天气温（摄氏度）的变量。

x = [32, 21, 24, 35, 10, 11, 22, 21, 27, 2]
y = [2000, 1000, 1500, 2500, 500, 900, 1100, 1500, 1800, 250]
df = pd.DataFrame([x,y]).T
df.columns = ['X', 'Y']
df
X Y
0 32 2000
1 21 1000
2 24 1500
3 35 2500
4 10 500
5 11 900
6 22 1100
7 21 1500
8 27 1800
9 2 250

我们可以通过散点图来观察 \(y_i\)（冰淇淋销售额（美元($’s)）和 \(x_i\)（摄氏度）之间的关系。

ax = df.plot(
    x='X', 
    y='Y', 
    kind='scatter', 
    ylabel=r'冰淇淋销售额（$\'s）', 
    xlabel='摄氏度'
)
_images/279b29dd784e43cef54622443b93d406a594a19a450a92702518347d4bb27aec.png

Fig. 44.1 散点图#

如你所见，数据表明在较热的日子里通常会卖出较多的冰淇淋。

为了构建数据的线性模型，我们需要选择代表“最佳”拟合线的 \(\alpha\) 和 \(\beta\) 值，使得

\[ \hat{y_i} = \hat{\alpha} + \hat{\beta} x_i \]

让我们从 \(\alpha = 5\) 和 \(\beta = 10\) 开始

α = 5
β = 10
df['Y_hat'] = α + β * df['X']

fig, ax = plt.subplots()
ax = df.plot(x='X',y='Y', kind='scatter', ax=ax)
ax = df.plot(x='X',y='Y_hat', kind='line', ax=ax, label=r'$\hat Y$')
plt.show()
_images/468063cda355bbddd56f1b13ddc1a67f28c47f78cf8cc8dfb85b8029108263a7.png

Fig. 44.2 带有拟合线的散点图#

我们可以看到这个模型没法很好地估计两者的关系。

我们可以继续猜测，并通过调整参数迭代出一条 “最佳 ”拟合线。

β = 100
df['Y_hat'] = α + β * df['X']

fig, ax = plt.subplots()
ax = df.plot(x='X',y='Y', kind='scatter', ax=ax)
ax = df.plot(x='X',y='Y_hat', kind='line', ax=ax, label=r'$\hat Y$')
plt.show()
_images/24a9ea5ecff2e64e5e4b2990e38578a7a2b9ec6050f665a73b5673d0e8da99cf.png

Fig. 44.3 带拟合线的散点图 #2#

β = 65
df['Y_hat'] = α + β * df['X']

fig, ax = plt.subplots()
ax = df.plot(x='X',y='Y', kind='scatter', ax=ax)
yax = df.plot(x='X',y='Y_hat', kind='line', ax=ax, color='g', label=r'$\hat Y$')
plt.show()
_images/f63b5c3ece86f59b82c8b04503b9cfdb737e92d07c925a576570bbcc7004eef1.png

Fig. 44.4 带拟合线的散点图 #3#

与其不断猜测参数值，我们可以把这个问题转化为一个优化问题，用数学方法来求解最优的参数。

为此，我们先来定义一个重要的概念：残差（residual）。残差 \(\epsilon_i\) 是实际观测值 \(y_i\) 与模型预测值 \(\hat{y}_i\) 之间的差异。

\[\begin{split} \begin{aligned} \hat{e}_i &= y_i - \hat{y}_i \\ &= y_i - \hat{\alpha} - \hat{\beta} x_i \end{aligned} \end{split}\]
df['error'] = df['Y_hat'] - df['Y']

df
X Y Y_hat error
0 32 2000 2085 85
1 21 1000 1370 370
2 24 1500 1565 65
3 35 2500 2280 -220
4 10 500 655 155
5 11 900 720 -180
6 22 1100 1435 335
7 21 1500 1370 -130
8 27 1800 1760 -40
9 2 250 135 -115 
fig, ax = plt.subplots()
ax = df.plot(x='X',y='Y', kind='scatter', ax=ax)
yax = df.plot(x='X',y='Y_hat', kind='line', ax=ax, color='g', label=r'$\hat Y$')
plt.vlines(df['X'], df['Y_hat'], df['Y'], color='r')
plt.show()
_images/9b40c1c9796b8cd275ca8bcee90a1dc046c31ef45605f0cd320d90a37b9122f0.png

Fig. 44.5 残差图#

普通最小二乘方法 (OLS) 通过 最小化 残差平方和 (SSR) 来选择 \(\alpha\) 和 \(\beta\) 的值。

\[ \min_{\alpha,\beta} \sum_{i=1}^{N}{\hat{e}_i^2} = \min_{\alpha,\beta} \sum_{i=1}^{N}{(y_i - \alpha - \beta x_i)^2} \]

我们称之为成本函数

\[ C = \sum_{i=1}^{N}{(y_i - \alpha - \beta x_i)^2} \]

我们希望通过调整参数 \(\alpha\) 和 \(\beta\) 来最小化这个成本函数。

44.1. 残差相对于 \(\alpha\) 和 \(\beta\) 是如何变化的#

首先，我们看看总残差相对于 \(\beta\) 的变化（保持截距 \(\alpha\) 不变）

我们从下一节可以知道 \(\alpha\) 和 \(\beta\) 的最优值是：

β_optimal = 64.38
α_optimal = -14.72

我们可以计算一系列 \(\beta\) 值的残差

errors = {}
for β in np.arange(20,100,0.5):
    errors[β] = abs((α_optimal + β * df['X']) - df['Y']).sum()

绘制残差图

ax = pd.Series(errors).plot(xlabel='β', ylabel='残差')
plt.axvline(β_optimal, color='r');
_images/233a641767a74e980b7239f1c74007309d7db6a2ef1ca9855e667efdeb374112.png

Fig. 44.6 绘制残差图#

现在我们改变 \(\alpha\) （保持 \(\beta\) 不变）

errors = {}
for α in np.arange(-500,500,5):
    errors[α] = abs((α + β_optimal * df['X']) - df['Y']).sum()

绘制残差图

ax = pd.Series(errors).plot(xlabel='α', ylabel='残差')
plt.axvline(α_optimal, color='r');
_images/edd915afb3bf4d63fe48fe30a95a3935a7df289c89b10f6a0d385c282d367870.png

Fig. 44.7 绘制残差图 (2)#

44.2. 计算最优值#

现在让我们使用微积分来解决优化问题，并计算出 \(\alpha\) 和 \(\beta\) 的最优值，以找到普通最小二乘（OLS）解。

首先对 \(\alpha\) 取偏导

\[ \frac{\partial C}{\partial \alpha}[\sum_{i=1}^{N}{(y_i - \alpha - \beta x_i)^2}] \]

并将其设为 \(0\)

\[ 0 = \sum_{i=1}^{N}{-2(y_i - \alpha - \beta x_i)} \]

我们可以通过两边除以 \(-2\) 来移除求和中的常数 \(-2\)

\[ 0 = \sum_{i=1}^{N}{(y_i - \alpha - \beta x_i)} \]

现在我们可以将这个方程分解成几个组成部分

\[ 0 = \sum_{i=1}^{N}{y_i} - \sum_{i=1}^{N}{\alpha} - \beta \sum_{i=1}^{N}{x_i} \]

中间项是常数 \(\alpha\) 从 \(i=1,...N\) 进行直接相加得到的

\[ 0 = \sum_{i=1}^{N}{y_i} - N*\alpha - \beta \sum_{i=1}^{N}{x_i} \]

重新排列各项可得

\[ \alpha = \frac{\sum_{i=1}^{N}{y_i} - \beta \sum_{i=1}^{N}{x_i}}{N} \]

我们可以发现分解成两个分数后，它们分别是均值 \(\bar{y_i}\) 和 \(\bar{x_i}\)

(44.1)#\[ \alpha = \bar{y_i} - \beta\bar{x_i} \]

回到成本函数 \(C\) ，现在我们对 \(\beta\) 取偏导

\[ \frac{\partial C}{\partial \beta}[\sum_{i=1}^{N}{(y_i - \alpha - \beta x_i)^2}] \]

并将其设为 \(0\)

\[ 0 = \sum_{i=1}^{N}{-2 x_i (y_i - \alpha - \beta x_i)} \]

我们可以通过将两边除以 \(-2\) ，再次将常数从求和中取出

\[ 0 = \sum_{i=1}^{N}{x_i (y_i - \alpha - \beta x_i)} \]

进一步化简

\[ 0 = \sum_{i=1}^{N}{(x_i y_i - \alpha x_i - \beta x_i^2)} \]

现在代入 \(\alpha\)

\[ 0 = \sum_{i=1}^{N}{(x_i y_i - (\bar{y_i} - \beta \bar{x_i}) x_i - \beta x_i^2)} \]

重新排列各项

\[ 0 = \sum_{i=1}^{N}{(x_i y_i - \bar{y_i} x_i - \beta \bar{x_i} x_i - \beta x_i^2)} \]

这个方程可以被分成两个求和

\[ 0 = \sum_{i=1}^{N}(x_i y_i - \bar{y_i} x_i) + \beta \sum_{i=1}^{N}(\bar{x_i} x_i - x_i^2) \]

得到 \(\beta\) 的解

(44.2)#\[ \beta = \frac{\sum_{i=1}^{N}(x_i y_i - \bar{y_i} x_i)}{\sum_{i=1}^{N}(x_i^2 - \bar{x_i} x_i)} \]

我们现在可以使用(44.1) 和 (44.2) 来计算\(\alpha\)和\(\beta\)的最优值

计算\(\beta\)

df = df[['X','Y']].copy()  # 原始数据

# 计算样本均值
x_bar = df['X'].mean()
y_bar = df['Y'].mean()

现在我们用10个观察值进行计算，然后把分子和分母分别求和

# 计算求和项
df['num'] = df['X'] * df['Y'] - y_bar * df['X']
df['den'] = pow(df['X'],2) - x_bar * df['X']
β = df['num'].sum() / df['den'].sum()
print(β)

64.37665782493369

计算\(\alpha\)

α = y_bar - β * x_bar
print(α)

-14.72148541114052

现在我们可以绘制OLS解

df['Y_hat'] = α + β * df['X']
df['error'] = df['Y_hat'] - df['Y']

fig, ax = plt.subplots()
ax = df.plot(x='X',y='Y', kind='scatter', ax=ax)
yax = df.plot(x='X',y='Y_hat', kind='line', ax=ax, color='g', label=r'$\hat Y$')
plt.vlines(df['X'], df['Y_hat'], df['Y'], color='r');
_images/b354f127d932a904b93d5c2bd44d6524ced7770ad67006ed938a340d65f0f721.png

Fig. 44.8 OLS最佳拟合线#

Exercise 44.1

现在你已经知道了使用OLS解决简单线性回归模型的方程，你可以开始运行自己的回归来构建\(y\)和\(x\)之间的模型了。

让我们考虑两个经济变量，人均GDP和预期寿命。

  1. 你认为它们之间的关系会是怎样的？

  2. 从our world in data中搜集一些数据

  3. 使用pandas导入csv格式的数据，并绘制几个不同国家的图表

  4. 使用(44.1) 和 (44.2)计算\(\alpha\)和\(\beta\)的最优值

  5. 使用OLS绘制最佳拟合线

  6. 解释系数并用一句话总结人均GDP和预期寿命之间的关系

Exercise 44.2

最小化平方和并不是生成最佳拟合线的 唯一 方法。

举个例子，我们还可以考虑最小化 绝对值 之和，这样可以减少对异常值的权重。

使用最小绝对值法求解 \(\alpha\) 和 \(\beta\)