求解线性回归系数
已知 n n n个观测值集合 { ( x i , y i ) , i = 1 , 2 , . . . , n } \{(x_i, y_i), i=1,2,...,n\} {(xi,yi),i=1,2,...,n}, 求回归系数 a a a,使得预测值 y ^ i = x i a \hat{y}_i={x_ia} y^i=xia与真实值 y i y_i yi的偏差平方和最小,即找目标函数 s = ∑ ( y i − y ^ ) 2 s=\sum(y_i - \hat{y})^2 s=∑(yi−y^)2的最小值。
- 当为一元线性回归,则 y i = x i 0 ∗ a 0 + x i 1 ∗ a 1 y_i = x_{i0}*a_0 + x_{i1}*a_1 yi=xi0∗a0+xi1∗a1,这里 x 0 i x_{0i} x0i恒等于1,那么 a 0 a_0 a0可看作为偏移量常数(截距);
- 当为多元( m m m元)线性回归时, x i , a i x_i,a_i xi,ai为向量,令 x i = ( x i 0 , x i 1 , x i 2 , . . . , x i m ) T \boldsymbol{x_i}= (x_{i0}, x_{i1}, x_{i2},...,x_{im})^T xi=(xi0,xi1,xi2,...,xim)T, a = ( a 0 , a 1 , a 2 , . . . , a m ) T \boldsymbol{a} = (a_{0}, a_{1}, a_{2},...,a_{m})^T a=(a0,a1,a2,...,am)T, 则 y i = ∑ j = 0 m x i j a i = x i T a {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
=
1
n
(
y
i
−
x
i
T
a
)
2
=
(
y
−
X
T
a
)
T
(
y
−
X
T
a
)
,
\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 = ( x 1 , x 2 , . . . , x n ) T \boldsymbol{X}=(\boldsymbol{x_{1},x_{2},...,x_{n}})^T X=(x1,x2,...,xn)T, y = ( y 1 , y 2 , . . . , y n ) T \boldsymbol{y}=(y_1,y_2,...,y_n)^T y=(y1,y2,...,yn)T.
求
s
\boldsymbol{s}
s的最小值,则可对目标函数
s
\boldsymbol{s}
s求导,令
u
=
y
−
X
T
a
\boldsymbol{u} = \boldsymbol{y} - \boldsymbol{X}^T\boldsymbol{a}
u=y−XTa:
s
′
(
a
)
=
(
u
T
u
)
′
=
u
T
u
′
+
u
T
u
′
=
2
u
T
(
−
X
T
)
=
−
2
(
X
u
)
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.
【标量对向量求导: ( u T v ) ′ = u T v ′ + v T u ′ (u^Tv)'=u^Tv'+v^Tu' (uTv)′=uTv′+vTu′】(u(x): nx1, v(x):nx1)
令 s ′ ( a ) = 0 \boldsymbol{s}'(\boldsymbol{a})=0 s′(a)=0,即 ( X ( y − X T a ) ) T = 0 (\boldsymbol{X} ( \boldsymbol{y} - \boldsymbol{X}^T\boldsymbol{a}))^T=0 (X(y−XTa))T=0,解得 a ^ = ( X X T ) − 1 X y \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^才有解。
上述方法求解回归系数是一般最小二乘法( o r d i n a r y l e a s t q u a r e s ordinary\ least\ quares ordinary least quares, OLS)
python实现
python中numpy中包含线性代数模块(linalg, linear algebra)可用于求解 a x = 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!') # 行列式为0 return None else: # 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 plt x_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 data regress_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)
附:相关系数计算公式:
C o r r ( X , Y ) = C o v ( X , Y ) V a r ( X ) V a r ( Y ) , Corr(X,Y) = \frac{Cov(X,Y)}{\sqrt{Var(X)}\sqrt{Var(Y)}}, Corr(X,Y)=Var(X)Var(Y)Cov(X,Y),
其中,- 协方差 C o v ( X , Y ) = E ( X Y ) − E ( X ) E ( Y ) Cov(X,Y)=E(XY) - E(X)E(Y) Cov(X,Y)=E(XY)−E(X)E(Y)
- V a r ( X ) , V a r ( Y ) Var(X), Var(Y) Var(X),Var(Y)分别为 X , Y X,Y X,Y的方差
- E ( X ) , E ( Y ) , E ( X Y ) E(X), E(Y), E(XY) E(X),E(Y),E(XY)分别为对用期望
- 协方差>0,则X与Y正相关;协方差<0,则负相关;协方差=0,则独立/不相关;同样相关系数与协方差同符号(同正负零),相关系数反应 X , Y X,Y X,Y的相关程度,其取值范围是 − 1 < C o r r ( X , Y ) < 1 -1<Corr(X, Y) < 1 −1<Corr(X,Y)<1, 即0<|Corr(X, Y)|<1,|Corr(X, Y)|的值越接近1,相关程度越高,反之,相关程度越低。