求解线性回归系数

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

• 当为一元线性回归，则$y_i=x_{i0}\ast a_0+x_{i1}\ast a_1$，这里$x_{0i}$恒等于1，那么$a_0$可看作为偏移量常数(截距)；
• 当为多元($m$元)线性回归时，$x_i,a_i$为向量，令${\mathbf{x}}_{\mathbf{i}}=(x_{i0},x_{i1},x_{i2},...,x_{im}{)}^T$, $\mathbf{a}=(a_0,a_1,a_2,...,a_m{)}^T$, 则$y_i={\sum }_{j=0}^m{x}_{ij}{a}_i={{\mathbf{x}}_{\mathbf{i}}}^T\mathbf{a}$.

因此，目标函数$\mathbf{s}$可用矩阵形式表示：
$$\mathbf{s}(\mathbf{a})=\sum _{i=1}^n(y_i-{{\mathbf{x}}_{\mathbf{i}}}^T\mathbf{a}{)}^2=(\mathbf{y}-{\mathbf{X}}^T\mathbf{a}{)}^T(\mathbf{y}-{\mathbf{X}}^T\mathbf{a}),$$

其中，$\mathbf{X}=({\mathbf{x}}_1,{\mathbf{x}}_2,...,{\mathbf{x}}_{\mathbf{n}}{)}^T$, $\mathbf{y}=(y_1,y_2,...,y_n{)}^T$.

求$\mathbf{s}$的最小值，则可对目标函数$\mathbf{s}$求导，令$\mathbf{u}=\mathbf{y}-{\mathbf{X}}^T\mathbf{a}$：
$${\mathbf{s}}^{\prime }(\mathbf{a})=({\mathbf{u}}^T\mathbf{u}{)}^{\prime }={\mathbf{u}}^T{\mathbf{u}}^{\prime }+{\mathbf{u}}^T{\mathbf{u}}^{\prime }=2{\mathbf{u}}^T(-{\mathbf{X}}^T)=-2(\mathbf{X}\mathbf{u}{)}^T.$$

【标量对向量求导：$(u^{T}v{)}^{\prime }=u^T{v}^{\prime }+v^T{u}^{\prime }$】(u(x): nx1, v(x):nx1)

令${\mathbf{s}}^{\prime }(\mathbf{a})=0$，即$(\mathbf{X}(\mathbf{y}-{\mathbf{X}}^T\mathbf{a}){)}^T=0$，解得 $\hat{\mathbf{a}}=({\mathbf{X}\mathbf{X}}^T{)}^{-1}\mathbf{X}\mathbf{y}$。 注意到，$\hat{\mathbf{a}}$的解中包含矩阵的逆，也就是说，只有当${\mathbf{X}}^{-1}$存在时，$\hat{\mathbf{a}}$才有解。

上述方法求解回归系数是一般最小二乘法（$ordinaryleastquares$, OLS）

python实现

python中numpy中包含线性代数模块(linalg, linear algebra)可用于求解$\mathbf{a}\mathbf{x}=\mathbf{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)
  

  
  附：相关系数计算公式:
  $$Corr(X,Y)=\frac{Cov(X,Y)}{\sqrt{Var(X)}\sqrt{Var(Y)}},$$
  其中，

  • 协方差 $Cov(X,Y)=E(XY)-E(X)E(Y)$
  • $Var(X),Var(Y)$分别为$X,Y$的方差
  • $E(X),E(Y),E(XY)$分别为对用期望
  • 协方差>0，则X与Y正相关；协方差<0，则负相关；协方差=0，则独立/不相关；同样相关系数与协方差同符号（同正负零），相关系数反应$X,Y$的相关程度，其取值范围是$-1<Corr(X,Y)<1$, 即0<|Corr(X, Y)|<1，|Corr(X, Y)|的值越接近1，相关程度越高，反之，相关程度越低。