多元線性回歸
多元線性回歸模型
實際中有很多問題是一個因變量與多個自變量成線性相關，我們可以用一個多元線性回歸方程來表示。
為了方便計算，我們將上式寫成矩陣形式:
Y = XW
- 假設自變量維度為N
- W為自變量的系數，下標0 - N
- X為自變量向量或矩陣，X維度為N,為了能和W0對應，X需要在第一行插入一個全是1的列。
- Y為因變量
那么問題就轉變成，已知樣本X矩陣以及對應的因變量Y的值，求出滿足方程的W，一般不存在一個W是整個樣本都能滿足方程，畢竟現實中的樣本有很多噪聲。最一般的求解W的方式是最小二乘法。
最小二乘法
我們希望求出的W是最接近線性方程的解的，最接近我們定義為殘差平方和最小，殘差的公式和殘差平方和的公式如下：
上面的公式用最小殘差平方和的方式導出的，還有一種思路用最大似然的方式也能推導出和這個一樣的公式，首先對模型進行一些假設：
- 誤差等方差不相干假設，即每個樣本的誤差期望為0，每個樣本的誤差方差都為相同值假設為σ
- 誤差密度函數為正態分布 e ~ N(0, σ^2)
簡單推導如下:
由此利用最大似然原理導出了和最小二乘一樣的公式。
最小二乘法求解
二次函數是個凸函數，極值點就是最小點。只需要求導數=0解出W即可。
模擬數據
我們這里用R語言模擬實踐一下，由於我們使用的矩陣運算，這個公式一元和多元都是兼容的，我們為了可視化方便一點，我們就用R語言自帶的women數據做一元線性回歸，和多元線性回歸的方式基本一樣。
women數據如下
> women
height weight
1 58 115
2 59 117
3 60 120
4 61 123
5 62 126
6 63 129
7 64 132
8 65 135
9 66 139
10 67 142
11 68 146
12 69 150
13 70 154
14 71 159
15 72 164
體重和身高具有線性關系，我們做一個散點圖可以看出來：
我們用最小二乘推導出來的公式計算w如下
X <- cbind(rep(1, nrow(women)), women$height)
X.T <- t(X)
Y <- women$weight
w <- solve(X.T %*% X) %*% X.T %*% Y
> w
[,1]
[1,] -87.51667
[2,] 3.45000
> lm.result <- lm(women$weight~women$height)
> lm.result
Call:
lm(formula = women$weight ~ women$height)
Coefficients:
(Intercept) women$height
-87.52 3.45
上面的R代碼w使我們利用公式計算出來的，下邊是R語言集成的線性回歸函數擬合出來的，可以看出我們的計算結果是正確的，lm的只是小數點取了兩位而已，將回歸出來的函數畫到圖中看下回歸的效果。
畫圖對應的R代碼如下，用R的感覺.....太飄逸了。
> png(file="chart2.png")
> plot(women$height, women$weight)
> lines(women$height, X %*% w)
> dev.off()
梯度下降法
除了用正規方程方式求解W，也可以用最常見的梯度下降法求得W，因為最小二乘是個凸函數，所以這里找到的極小點就是最小點。下面這段代碼用R寫還是非常容易的，但是剛開始step步長參數調的太大了，導致一直不收斂，我還
以為是程序錯誤，后來怎么看也沒寫錯，就把參數調了個很小值，結果就收斂了。step的這個取值其實應該是變化的，先大后下比較科學，我這個調的很小，需要接近500萬次才能收斂。
- 初始化W 為全0向量，也可以隨機一個向量
- 設置最大迭代次數，本例為了收斂設置了一個很大的數
- 設置步長step，小了收斂很慢，大了不收斂.......
- 求損失函數的梯度
- W(k+1) 為 W(k) + 損失函數負梯度 * 步長step
- 循環，直到梯度接近0
X <- cbind(rep(1, nrow(women)), women$height)
Y <- women$weight
maxIterNum <- 5000000;
step <- 0.00003;
W <- rep(0, ncol(X))
for (i in 1:maxIterNum){
grad <- t(X) %*% (X %*% W - Y);
if (sqrt(as.numeric(t(grad) %*% grad)) < 1e-3){
print(sprintf('iter times=%d', i));
break;
}
W <- W - grad * step;
}
print(W);
輸出
[1] "iter times=4376771"
print(W);
[,1]
[1,] -87.501509
[2,] 3.449768
歸一化
上面的批量梯度下降為什么收斂如此之慢呢？原因很簡單，沒有做歸一化，做了歸一化，收斂速度快了非常非常多！！！！
正確代碼如下：
XScale = scale(women$height)
Ux = attr(XScale, "scaled:center")
Dx = attr(XScale, "scaled:scale")
YScale = scale(women$weight)
Uy = attr(YScale, "scaled:center")
Dy = attr(YScale, "scaled:scale")
X <- cbind(rep(1, nrow(women)), as.matrix(XScale))
Y <- as.matrix(YScale)
maxIterNum <- 5000;
step <- 0.001;
W <- rep(0, ncol(X))
for (i in 1:maxIterNum){
grad <- t(X) %*% (X %*% W - Y);
if (sqrt(as.numeric(t(grad) %*% grad)) < 1e-6){
print(sprintf('iter times=%d', i));
break;
}
W <- W - grad * step;
}
print(W);
W0 = W[1]
Wn = W[2:length(W)]
Wn = Dy * Wn / Dx
W0 = Uy + Dy * W0 - Dy * Ux / Dx
W = c(W0, Wn)
print(W);
輸出
[1] "iter times=1168"
print(W);
-88.53154 3.45000
更多精彩文章 http://h2cloud.org/