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()

一元线性回归的Python实现