线性回归模型和最小二乘法

线性回归模型和最小二乘法

最小二乘法极小化残差的平方和，该准则度量平均拟合偏离。
将残差平方和写成如下形式
$RSS(\theta)=(y-X\beta)^T(y-X\beta)$
这是$p+1$个参数的二次函数。
关于$\beta$微分，得到
${\partial RSS \over \partial \beta} = -2X^T(y-X^T\beta)$
假定$X$是列满秩的，从而$X^TX$是正定的。
令$X^T(y-X^T\beta)=0$
得到唯一解：
$\hat \beta = (X^TX)^{-1}X^Ty$
在训练输入上的拟合值为
$\hat y = X\hat \beta = X(X^TX)^{-1}y$
从几何上来看，这个拟合向量是y向X的列空间的投影，矩阵$P = X(X^TX)X^T$也称为投影矩阵，同时也是一个幂等矩阵，即满足
$P^n = P$
为了确定$\hat \beta$的性质，假定观测$y_i$是不相关的，具有常数方差$\sigma ^2$，而$x_i$是固定的，非随机的。
即y满足$y=x\beta+\epsilon,\epsilon$为观测误差服从正态分布$N(0,\sigma ^2)$
由此可以得到$\hat \beta$的方差-协方差矩阵
$Var(\hat \beta) = Var((X^TX)^{-1}X^Ty) = (X^TX)^{-1}\sigma ^2$
现在来求$\sigma ^2$的估计
$\sigma ^2$等于误差e的平方期望值，我们想使用残差的平方和来估计方差$\sigma ^2$,
残差向量为
$\hat e = y - \hat y = y - X(X^TX)^{-1}X^Ty$
将投影矩阵用P表示，P有如下性质
$P = P^T = P^n$
$(I-P) = (I-P)^T = (I-P)^n$
从几何上来理解，将一个向量在一个线性空间内投影多次等同于投影一次。
残差的平方和为
$||y-Py||^2 = ||(I-P)y||^2=y^T(I-P)(I-P)y=y^T(I-P)y$
这里用到了一个定理：
x是n维随机变量，具有均值$u$和协方差矩阵$\Sigma$，A是一个固定矩阵，那么
$E(x^TAx) = tr(A\Sigma)+u^TAu$
那么我们可以得到：
$E(y^T(I-P)y) = \sigma ^2tr(I-P) + E(y^T)*(I-P)*E(y)$
因为$E(y) = X\beta$，所以$(I-P)*E(y)=X\beta-X(X^TX)^{-1}X^TX\beta=0$
此外$tr(P) = tr(X(X^TX)^{-1}X^T)$,
利用迹的交换律，得到$tr(P) = tr(X^TX(X^TX)^{-1})=tr(I_{(p+1)*（p+1)})=p+1$
$tr(I-P)=tr(I_{N*N})-tr(P)=n-p-1$,
所以$E(||y-Py||^2)=E(y^T(I-P)y)=\sigma ^2(n-p)$
所以可以得到方差的无偏估计
$\hat \sigma ^2 = ||y-Py||^2 / (n-p+1)$
然后可以得到以下性质（证明略）
$\hat \beta$服从高斯分布$N(\beta,(X^TX)^{-1}\sigma ^2)$
${（N-p-1)\hat \sigma ^2\over \sigma ^2}$服从自由度为N-p-1的卡方分布
$\hat \beta$和$\hat \sigma ^2$是统计独立的

利用这些分布性质，可以形成参数$\beta_j$的假设检验和置信区间
为检验特定系数$\beta_j=0$的假设，形成Z得分
$z_j = {\hat \beta \over \sigma \sqrt{v_j}}$
其中$v_j$是$(X^TX)^{-1}$的第j个对角线元素。在原假设$\beta_j=0$下，$z_j$服从分布$t_{n-p-1}$（具有自由度为$N-p-1$的t分布），
因此绝对值大的$z_j$将导致拒绝该假设
当想要检验一组变量能否从一个模型里排除的时候，可以使用F统计量
$F={(RSS_0-RSS_1)/(p_1-p_0) \over RSS1/(N-p1-1)}$
其中$RSS_1$是具有$p_1+1$个参数的较大模型的残差，$RSS_0$是具有$p_0+1$的变量的残差和。
F统计量具有分布$F_{p_1-p_0,N-p_1-1}$,可以证明当删除一个系数时，F统计量相当于Z得分的平方。
另外也可以通过类似的方法得到$\beta$的置信区间。

高斯-马尔科夫定理

高斯马尔科夫定理：在所有的线性无偏估计中，参数$\beta$的最小二乘方估计$\hat \beta$具有最小方差。
也就是说假设输入为$\alpha$，如果有$\alpha^T\beta$的其他无偏估计$cy$，即
$E(cy)=\alpha^T\beta$
那么$Var(\alpha^T\hat \beta) \le Var(cy)$
现在简单的证明一下高斯-马尔科夫定理。
$c$可以写成$\alpha^T(X^TX)^{-1}X^T+d$的形式
然后根据
$E(cy)=cE(y)=(\alpha^T(X^TX)^{-1}X^T+d)X\beta$
$=\alpha\beta+dX\beta=\alpha^T\beta$
对于任意的$\beta$，有
$dx\beta=0$，所以得到
$dx=0$
首先求出最小二乘法估计的方差
$Var(\alpha^T\hat \beta)=\alpha^TVar(\hat \beta)\alpha$
$=\sigma ^2 \alpha^T(X^TX)^{-1}\alpha$
然后是其余估计的方差
$Var(cy)=cVar(y)c^T=\sigma ^2(\alpha^T(X^TX)^{-1}\alpha+dd^T)$
而dd^T是一个非负整数，所以得到
$Var(\alpha^T\hat \beta) \le Var(cy)$

代码实现

from numpy import *


# 预处理数据
def loadData(filename):
    dataSet = []
    labels = []
    fr = open(filename)
    for line in fr.readlines():
        curLine = line.strip().split('\t')
        fltLine = list(map(float, curLine))
        dataSet.append(fltLine[:-1])
        labels.append(fltLine[-1])
    return dataSet, labels



# 普通最小二乘法
def ols(xArr, yArr):
    xMat = mat(xArr); yMat = mat(yArr).T
    xTx = xMat.T*xMat
    if linalg.det(xTx) == 0.0:
        print("xTx is not invertible")
        return
    return xTx.I * (xMat.T) * yMat