線性回歸模型和最小二乘法

最小二乘法極小化殘差的平方和，該準則度量平均拟合偏離。

将殘差平方和寫成如下形式

RSS(θ)=(y−Xβ)T(y−Xβ)

這是 p+1 個參數的二次函數。

關于 β 微分，得到

∂RSS∂β=−2XT(y−XTβ)

假定 X 是列滿秩的，進而XTX是正定的。

令 XT(y−XTβ)=0

得到唯一解：

β^=(XTX)−1XTy

在訓練輸入上的拟合值為

y^=Xβ^=X(XTX)−1y

從幾何上來看，這個拟合向量是y向X的列空間的投影，矩陣 P=X(XTX)XT 也稱為投影矩陣，同時也是一個幂等矩陣，即滿足

Pn=P

為了确定 β^ 的性質，假定觀測 yi 是不相關的，具有常數方差 σ2 ，而 xi 是固定的，非随機的。

即y滿足 y=xβ+ϵ,ϵ 為觀測誤差服從正态分布 N(0,σ2)

由此可以得到 β^ 的方差-協方差矩陣

Var(β^)=Var((XTX)−1XTy)=(XTX)−1σ2

現在來求 σ2 的估計

σ2 等于誤差e的平方期望值，我們想使用殘差的平方和來估計方差 σ2 ,

殘差向量為

e^=y−y^=y−X(XTX)−1XTy

将投影矩陣用P表示，P有如下性質

P=PT=Pn

(I−P)=(I−P)T=(I−P)n

從幾何上來了解，将一個向量在一個線性空間内投影多次等同于投影一次。

殘差的平方和為

||y−Py||2=||(I−P)y||2=yT(I−P)(I−P)y=yT(I−P)y

這裡用到了一個定理：

x是n維随機變量，具有均值 u 和協方差矩陣Σ，A是一個固定矩陣，那麼

E(xTAx)=tr(AΣ)+uTAu

那麼我們可以得到：

E(yT(I−P)y)=σ2tr(I−P)+E(yT)∗(I−P)∗E(y)

因為 E(y)=Xβ ，是以 (I−P)∗E(y)=Xβ−X(XTX)−1XTXβ=0

此外 tr(P)=tr(X(XTX)−1XT) ,

利用迹的交換律，得到 tr(P)=tr(XTX(XTX)−1)=tr(I(p+1)∗（p+1))=p+1

tr(I−P)=tr(IN∗N)−tr(P)=n−p−1 ,

是以 E(||y−Py||2)=E(yT(I−P)y)=σ2(n−p)

是以可以得到方差的無偏估計

σ^2=||y−Py||2/(n−p+1)

然後可以得到以下性質（證明略）

β^ 服從高斯分布 N(β,(XTX)−1σ2)

（N−p−1)σ^2σ2 服從自由度為N-p-1的卡方分布

β^ 和 σ^2 是統計獨立的

利用這些分布性質，可以形成參數 βj 的假設檢驗和置信區間

為檢驗特定系數 βj=0 的假設，形成Z得分

zj=β^σvj√

其中 vj 是 (XTX)−1 的第j個對角線元素。在原假設 βj=0 下， zj 服從分布 tn−p−1 （具有自由度為 N−p−1 的t分布），

是以絕對值大的 zj 将導緻拒絕該假設

當想要檢驗一組變量能否從一個模型裡排除的時候，可以使用F統計量

F=(RSS0−RSS1)/(p1−p0)RSS1/(N−p1−1)

其中 RSS1 是具有 p1+1 個參數的較大模型的殘差， RSS0 是具有 p0+1 的變量的殘差和。

F統計量具有分布 Fp1−p0,N−p1−1 ,可以證明當删除一個系數時，F統計量相當于Z得分的平方。

另外也可以通過類似的方法得到 β 的置信區間。

高斯-馬爾科夫定理

高斯馬爾科夫定理：在所有的線性無偏估計中，參數 β 的最小二乘方估計 β^ 具有最小方差。

也就是說假設輸入為 α ，如果有 αTβ 的其他無偏估計 cy ，即

E(cy)=αTβ

那麼 Var(αTβ^)≤Var(cy)

現在簡單的證明一下高斯-馬爾科夫定理。

c 可以寫成αT(XTX)−1XT+d的形式

然後根據

E(cy)=cE(y)=(αT(XTX)−1XT+d)Xβ

=\alpha\beta+dX\beta=\alpha^T\beta =αβ+dXβ=αTβ

對于任意的\beta β ，有

dx\beta=0 dxβ=0 ，是以得到

dx=0 dx=0

首先求出最小二乘法估計的方差

Var(\alpha^T\hat \beta)=\alpha^TVar(\hat \beta)\alpha Var(αTβ^)=αTVar(β^)α

=σ2αT(XTX)−1α

然後是其餘估計的方差

Var(cy)=cVar(y)cT=σ2(αT(XTX)−1α+ddT)

而dd^T是一個非負整數，是以得到

Var(αTβ^)≤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[:-])
        labels.append(fltLine[-])
    return dataSet, labels



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