线性回归系数求解及Python实现

求解线性回归系数

已知nnn个观测值集合{(xi,yi),i=1,2,...,n}\{(x_i, y_i), i=1,2,...,n\}{(xi,yi),i=1,2,...,n}, 求回归系数aaa，使得预测值y^i=xia\hat{y}_i={x_ia}y^i=xia与真实值yiy_iyi的偏差平方和最小，即找目标函数s=∑(yi−y^)2s=\sum(y_i - \hat{y})^2s=∑(yi−y^)2的最小值。

当为一元线性回归，则yi=xi0∗a0+xi1∗a1y_i = x_{i0}*a_0 + x_{i1}*a_1yi=xi0∗a0+xi1∗a1，这里x0ix_{0i}x0i恒等于1，那么a0a_0a0可看作为偏移量常数(截距)；
当为多元(mmm元)线性回归时，xi,aix_i,a_ixi,ai为向量，令xi=(xi0,xi1,xi2,...,xim)T\boldsymbol{x_i}= (x_{i0}, x_{i1}, x_{i2},...,x_{im})^Txi=(xi0,xi1,xi2,...,xim)T, a=(a0,a1,a2,...,am)T\boldsymbol{a} = (a_{0}, a_{1}, a_{2},...,a_{m})^Ta=(a0,a1,a2,...,am)T, 则yi=∑j=0mxijai=xiTa{y_i}=\sum_{j=0}^m{x_{ij}a_i} =\boldsymbol{x_i}^T\boldsymbol{a}yi=∑j=0mxijai=xiTa.

因此，目标函数s\boldsymbol{s}s可用矩阵形式表示：
s(a)=∑i=1n(yi−xiTa)2=(y−XTa)T(y−XTa),\boldsymbol{s}(\boldsymbol{a})=\sum_{i=1}^n{({y_i}-\boldsymbol{x_i}^T\boldsymbol{a})^2} =(\boldsymbol{y} - \boldsymbol{X}^T\boldsymbol{a})^T (\boldsymbol{y} - \boldsymbol{X}^T\boldsymbol{a}),s(a)=i=1∑n(yi−xiTa)2=(y−XTa)T(y−XTa),

其中，X=(x1,x2,...,xn)T\boldsymbol{X}=(\boldsymbol{x_{1},x_{2},...,x_{n}})^TX=(x1,x2,...,xn)T, y=(y1,y2,...,yn)T\boldsymbol{y}=(y_1,y_2,...,y_n)^Ty=(y1,y2,...,yn)T.

求s\boldsymbol{s}s的最小值，则可对目标函数s\boldsymbol{s}s求导，令u=y−XTa\boldsymbol{u} = \boldsymbol{y} - \boldsymbol{X}^T\boldsymbol{a}u=y−XTa：
s′(a)=(uTu)′=uTu′+uTu′=2uT(−XT)=−2(Xu)T.\boldsymbol{s}'(\boldsymbol{a}) = (\boldsymbol{u}^T\boldsymbol{u})' =\boldsymbol{u}^T\boldsymbol{u}'+\boldsymbol{u}^T\boldsymbol{u}'=2\boldsymbol{u}^T(-\boldsymbol{X}^T)=-2(\boldsymbol{X}\boldsymbol{u})^T.s′(a)=(uTu)′=uTu′+uTu′=2uT(−XT)=−2(Xu)T.

【标量对向量求导：(uTv)′=uTv′+vTu′(u^Tv)'=u^Tv'+v^Tu'(uTv)′=uTv′+vTu′】(u(x): nx1, v(x):nx1)

令s′(a)=0\boldsymbol{s}'(\boldsymbol{a})=0s′(a)=0，即(X(y−XTa))T=0(\boldsymbol{X} ( \boldsymbol{y} - \boldsymbol{X}^T\boldsymbol{a}))^T=0(X(y−XTa))T=0，解得 a^=(XXT)−1Xy\hat{\boldsymbol{a}}=(\boldsymbol{XX}^T)^{-1}\boldsymbol{X}\boldsymbol{y}a^=(XXT)−1Xy。注意到，a^\hat{\boldsymbol{a}}a^的解中包含矩阵的逆，也就是说，只有当X−1\boldsymbol{X}^{-1}X−1存在时，a^\hat{\boldsymbol{a}}a^才有解。

上述方法求解回归系数是一般最小二乘法（ordinaryleastquaresordinary\ least\ quaresordinary least quares, OLS）

python实现

python中numpy中包含线性代数模块(linalg, linear algebra)可用于求解ax=b\boldsymbol{ax=b}ax=b。

一、导入数据

from numpy import *
import pandas as pd
stu_score = '{mydata_path}/data.tsv'
dataset = pd.read_csv(stu_score, index_col=False)
dataset.head()
len(dataset)

二、数据转化为矩阵，并计算回归系数

def load_data(data_file):data_arr = []label_arr = []with open(data_file, 'r') as f:header = f.readline()for line in f:mydata = line.strip().split(',')mydata.insert(0, 1)  # 假定偏移量是常数c，第一列补1(y = ax + c)data_info = [float(i) for i in mydata]data_arr.append(data_info[:-1])label_arr.append([data_info[-1]])  # 最后一列为对应y值return mat(data_arr), mat(label_arr)def stand_regres(x_mat, y_mat):x_mat_T = x_mat.T * x_mat # 下面判断x_mat_T是否可逆if linalg.det(x_mat_T) == 0:  # 若行列式|x_mat_T|不为0，则A可逆print('This matrix is singular, cannot do inverse!')  # 行列式为0return Noneelse:# reg_coef = x_mat_T.I * (x_mat.T * y_mat)  # 可根据上面推到得到回归系数，reg_coef = linalg.solve(x_mat_T, x_mat.T * y_mat)  # 也可根据numpy中linalg模块中solve方法解ax + b = 0得到回归系数return reg_coef

三、数据可视化

import matplotlib.pyplot as pltx_mat, y_mat = load_data(stu_score)  # header
fig = plt.figure()
ax = fig.add_subplot(111)
ax.scatter(x_mat[:,1].flatten().A[0], y_mat[:,0].flatten().A[0])  # plot scatter of original dataregress_coef = stand_regres(x_mat, y_mat)
x_copy = x_mat.copy()
y_hat = x_copy * regress_coef
ax.plot(x_copy[:,1], y_hat)  # plot regression line
plt.show()

四、拟合直线的相关系数
```
corrcoef(y_hat.T, y_mat.T)
```
附：相关系数计算公式:
Corr(X,Y)=Cov(X,Y)Var(X)Var(Y),Corr(X,Y) = \frac{Cov(X,Y)}{\sqrt{Var(X)}\sqrt{Var(Y)}},Corr(X,Y)=Var(X)Var(Y)Cov(X,Y),
其中，
- 协方差 Cov(X,Y)=E(XY)−E(X)E(Y)Cov(X,Y)=E(XY) - E(X)E(Y)Cov(X,Y)=E(XY)−E(X)E(Y)
- Var(X),Var(Y)Var(X), Var(Y)Var(X),Var(Y)分别为X,YX,YX,Y的方差
- E(X),E(Y),E(XY)E(X), E(Y), E(XY)E(X),E(Y),E(XY)分别为对用期望
- 协方差>0，则X与Y正相关；协方差<0，则负相关；协方差=0，则独立/不相关；同样相关系数与协方差同符号（同正负零），相关系数反应X,YX,YX,Y的相关程度，其取值范围是−1<Corr(X,Y)<1-1<Corr(X, Y) < 1−1<Corr(X,Y)<1, 即0<|Corr(X, Y)|<1，|Corr(X, Y)|的值越接近1，相关程度越高，反之，相关程度越低。