一元线性回归的Python实现

目录

• 1 问题的提出
• 2 原理
  • 2.1 代价函数
  • 2.2 模型的评价
    • 2.2.1 皮尔逊相关系数
    • 2.2.2 决定系数
• 3 Python 实现
  • 3.1 不调sklearn库
  • 3.2 调 sklearn 库
• 4 梯度下降法
  • 4.1 原理
  • 4.2 Python实现
• 参考

1 问题的提出

对于给定的数据集 $D = \{(x_1,y_1),(x_2,y_2),\cdots,(x_m,y_m)\}$，线性回归 (linear regression) 试图学得一个线性模型以尽可能准确地预测是指输出标记.

2 原理

设给定的数据集 $D = \{(x_i,y_i)\}_{i=1}^m, \ x_i,y_i \in \mathcal{R}$. 对于离散属性，如果属性值间存在“序”(order)的关系，可通过连续化将其转化为连续值，例如二值属性“身高”的取值“高”“矮”可转化为 $\{1.0,0.0\}$，三值属性“高度”的取值“高”“中”“低”可转化为 $\{1.0,0.5,0.0\}$；若属性之间不存在有序关系，假定有 $k$ 个属性值，则通常转化为 $k$ 维向量，例如属性“瓜类”的取值“西瓜”“南瓜”“黄瓜”可转化为 $(0,0,1),(0,1,0),(1,0,0)$.

线性回归试图学得

$$f(x_i) = wx_i + b, \ s.t. f(x_i) \simeq y_i \tag{1}$$

2.1 代价函数

我们需要确定 $w$ 和 $b$ 的值使得 $f(x)$ 和 $y$ 之间的差别尽可能小. 因此我们引入均方误差，均方误差是回归问题中最常用的性能度量工具，我们可以试图让均方误差最小化，即

$$\begin{aligned} (w^*,b^*) &= \arg\min\limits_{(w,b)} \sum_{i=1}^m(f(x_i)-y_i)^2\\ &=\arg\min\limits_{(w,b)} \sum_{i=1}^m(y_i-wx_i-b)^2 \end{aligned} \tag{2}$$

$\arg\min\limits_{\theta}f(x;\theta)$ 意思就是找出一个 $\theta$ 使得 $f(x;\theta)$ 的值最小，即他的返回值是 $f(x;\theta)$ 的最小值所对应的 $\theta$ 的值.

均方误差有很好的几何意义，它对应了常用的“欧式距离”(Euclidean distance)，基于均方误差最小化的进行模型求解的方法称为“最小二乘法”(least square method). 在线性回归中，最小二乘法就是试图寻找一条直线，使得所有样本点到直线的欧氏距离之和最小.

在均方误差的基础上进一步构造代价函数

$$J(w,b) = \frac{1}{2m}\sum_{i=1}^m(f(x_i)-y_i)^2 = \frac{1}{2m}\sum_{i=1}^m(wx_i + b - y_i)^2$$

这里分母的 $2$ 是为了后续求导的方便

求解 $w$ 和 $b$ 使 $J(w,b)$ 最小化的过程，称为线性回归模型的最小二乘“参数估计”(parameter estimation). 我们可将 $J(w,b)$ 分别对 $w$ 和 $b$ 求导，得到

$$\begin{aligned} &\frac{\partial J(w,b)}{\partial w} = \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i)x_i\\ &\frac{\partial J(w,b)}{\partial b} = \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i) \end{aligned}$$

令上两式等于0，解得

$$w = \frac{\sum_{i=1}^my_i(x_i-\bar{x})}{\sum_{i=1}^mx_i^2-\frac{1}{m}(\sum_{i=1}^mx_i)^2}\\ b = \frac{1}{m}\sum_{i=1}^m(y_i-wx_i)$$

$J(w,b)$ 是关于 $w,b$ 的凸函数，根据凸函数的性质，其偏导为 $0$ 时就是 $w$ 和 $b$ 的最优解.

其中 $\bar{x} = \frac{1}{m}\sum_\limits{i=1}^mx_i$ 为 $x$ 的均值.

2.2 模型的评价

2.2.1 皮尔逊相关系数

使用相关系数衡量线性相关性的强弱，皮尔逊相关系数的公式如下：

$$r_{xy} = \frac{COV(X,Y)}{\sqrt{Var(X)Var(Y)}}$$

相关度越高，皮尔逊相关系数的值就趋于 1 或 -1 （趋于 1 表示它们呈正相关， 趋于 -1 表示它们呈负相关）；如果相关系数等于0，表明它们之间不存在线性相关关系.

2.2.2 决定系数

决定系数 $R^2$ 也称拟合优度，反应了 $y$ 的波动有多少百分比能被 $x$ 的波动所描述. 决定系数越接近 1 ，说明拟合程度越好.

总平方和

$$SST = \sum_{i=1}^n(y_i-\bar{y})^2$$

回归平方和

$$SSR = \sum_{i=1}^n(\hat{y}_i-\bar{y})^2$$

残差平方和

$$SSE = \sum_{i=1}^n(y_i-\hat{y}_i)^2$$

其中

$$SST = SSR + SSE\\ R^2 = \frac{SSR}{SST} = 1- \frac{SSE}{SST}$$

3 Python 实现

3.1 不调sklearn库

Step1. 调库

import numpy as np
from numpy import genfromtxt
import matplotlib.pyplot as plt

Step2. 数据导入并绘制散点图

data = genfromtxt("data.csv", delimiter = ",")
x = data[:, 0, np.newaxis]
y = data[:, 1, np.newaxis]
plt.scatter(x, y)
plt.show()

Step3. 求回归系数
根据先前的推导，已经知道

$$w = \frac{\sum_{i=1}^my_i(x_i-\bar{x})}{\sum_{i=1}^mx_i^2-\frac{1}{m}(\sum_{i=1}^mx_i)^2}\\ b = \frac{1}{m}\sum_{i=1}^m(y_i-wx_i)$$

m = len(x)
x_bar = np.mean(x)
w = (np.sum((x - x_bar)*y))/(np.sum(x**2)-(1/m)*(np.sum(x))**2)
b = np.mean(y-w*x)
print(w,b)

1.3224310227553517 7.991020982270779

Step4. 拟合图像

plt.plot(x, y, 'b.')
plt.plot(x, w*x+b, 'r')
plt.show()

Step5. 计算相关系数和决定系数

COVxy = np.cov(x.T,y.T)
r = COVxy[0,1]/(x.std()*y.std())
print(r)

0.78154393928063

y_hat = w*x+b
SSR = np.sum((y_hat-np.mean(y))**2)
SST = np.sum((y-np.mean(y))**2)
R2 = SSR/SST
print(R2)

0.5986557915386548

3.2 调 sklearn 库

建模：

from sklearn.linear_model import LinearRegression
model = LinearRegression()
model.fit(x, y)

LinearRegression()

拟合图像的得出

plt.plot(x, y, 'b.')
plt.plot(x, model.predict(x), 'r')
plt.show()

回归系数

w = model.intercept_
b = model.coef_
print("截距为 {0}, 回归系数为 {1}".format(w[0], b[0][0]))

截距为 7.991020982270385, 回归系数为 1.32243102275536

决定系数

model.score(x, y)

0.598655791538662

4 梯度下降法

4.1 原理

由于代价函数是凸函数，因此只有全局最小值，梯度下降法的原理是先定一个初始值，然后利用导数就像阶梯一样慢慢逼近全局最小值

图片出处：https://zhuanlan.zhihu.com/p/36564434

已知代价函数 $J(w,b)$，我们需要找一组 $w,b$ 使得 $J(w,b)$ 最小，给定一个算法:

$$\begin{aligned} &给定 \ w,b \ 的初始值 \ w_0,b_0\\ \\ &repeat until convergence \ \{\\ &\quad\quad temp0 = w - \alpha \frac{\partial}{\partial w}J(w,b)\\ &\quad\quad temp1 = b - \alpha \frac{\partial}{\partial b}J(w,b)\\ &\quad\quad w = temp0\\ &\} \end{aligned}$$

其中 $\alpha$ 为学习率，学习率不能太大也不能太小，可以多次尝试 $0.1,0.03,0.01,0.003,0.001,0.0003,\cdots$.

根据已知条件，有

$$J(w,b) = \frac{1}{2m}\sum_{i=1}^m(wx_i + b - y_i)^2\\ \frac{\partial J(w,b)}{\partial w} = \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i)x_i\\ \frac{\partial J(w,b)}{\partial b} = \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i)$$

于是

$$w = w - \alpha \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i)\\ b = b - \alpha \frac{1}{m}\sum_{i=1}^m(wx_i+b-y_i)x_i$$

4.2 Python实现

# 学习率learning rate
lr = 0.0001
# 截距初值
b = 0
# 斜率初值
w = 0
# 最大迭代次数
epochs = 50
# 最小二乘法
def compute_error(b, w, x, y):
  totalError = 0
  for i in range(0, len(x)):
    totalError += (y[i] - (w * x[i] + b)) ** 2
  return totalError / float(len(x)) / 2.0
 
def gradient_descent_runner(x, y, b, w, lr, epochs):
  # 计算总数据量
  m = float(len(x))
  # 循环epochs次
  for i in range(epochs):
    b_grad = 0
    w_grad = 0
    # 计算梯度的总和再求平均
    for j in range(0, len(x)):
      b_grad += -(1/m) * (y[j] - ((w * x[j]) + b))
      w_grad += -(1/m) * x[j] * (y[j] - ((w * x[j]) + b))
     
    # 更新 b 和 w
    b = b - (lr * b_grad)
    w = w - (lr * w_grad)
    
    # 每迭代5次，输出一次图像
    if i % 5 == 0:
      print('epochs:', i)
      plt.plot(x, y, 'b.')
      plt.plot(x, w*x + b, 'r')
      plt.show()
  return b,w

print('Starting b = {0}, w = {1}, error = {2}'.format(b, w, compute_error(b, w, x, y)))
print('Running')
b, w = gradient_descent_runner(x, y, b, w, lr, epochs)
print('After {0} iterations b = {1}, w = {2}, error = {3}'.format(epochs, b, w, compute_error(b,w,x,y)))
# 画图
plt.plot(x, y, 'b.')
plt.plot(x, w * x + b, 'r')
plt.show()

迭代50次后得到 $b = 0.03207192, w = 1.47886174$，最小二乘误差为 56.3244305.

参考

[1]. 周志华.《机器学习》.清华大学出版社
[2]. https://www.bilibili.com/video/BV1Rt411q7WJ?p=4&vd_source=08d2535b05740d396ec0fc720c52f36a