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^*) &= argminlimits_{(w,b)} sum_{i=1}^m(f(x_i)-y_i)^2\
&=argminlimits_{(w,b)} sum_{i=1}^m(y_i-wx_i-b)^2
end{aligned} tag{2}
]

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

一元线性回归的Python实现

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

一元线性回归的Python实现

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

一元线性回归的Python实现

回归系数

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 原理

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

图片出处：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 {\
&quadquad temp0 = w - alpha frac{partial}{partial w}J(w,b)\
&quadquad temp1 = b - alpha frac{partial}{partial b}J(w,b)\
&quadquad 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实现