一元线性回归的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()