【預測演算法】01. 線性回歸原理及R語言實現

01-27

回歸分析是統計學的核心演算法，是機器學習最基本演算法，也是數學建模最常用的演算法之一。

簡單來說，回歸分析就是用一個或多個自變數來預測因變數的方法，具體是通過多組自變數和因變數的樣本數據，擬合出最佳的函數關係。

本篇由前入深將線性回歸的原理講清楚，並用案例演示實際操作。

一、最小二乘法

設有 $n$ 組樣本點： $(x_i, y_i), quad i=1, cdots, n$

例1，現有10期的廣告費用與銷售額的數據：

先畫散點圖觀察一下：

cost<-c(30,40,40,50,60,70,70,70,80,90)nsale<-c(143.5,192.2,204.7,266,318.2,457,333.8,312.1,386.4,503.9)ndat<-as.data.frame(cbind(cost,sale))nplot(dat)n

可見，這些散點大致在一條直線上，一元線性回歸就是尋找一條直線，使得與這些散點擬合程度最好（越接近直線越好）。

比如畫這樣一條直線，方程可寫為： $y=beta_0 + beta_1 x$ （線性模型），其中 $beta_0, beta_1$ 是待定係數，目標是選取與樣本點最接近的直線對應的 $beta_0, beta_1$ .

那麼，怎麼刻畫這種「最接近」？

$hat{y}_i = beta_0 + beta_1 x_i$ 是與橫軸 $x_i$ 對應的直線上的點的縱坐標（稱為線性模型預測值），它與樣本點 $x_i$ 對應的真實值 $y_i$ 之差，就是預測誤差（紅線長度）：

$varepsilon_i = |y_i - hat{y}_i|, quad i = 1, cdots, n$

適合描述散點到直線的「接近程度」。

但絕對值不容易計算，改用：

$varepsilon_i^2 = (y_i - hat{y}_i)^2, quad i = 1, cdots, n$

我們需要讓所有散點總體上最接近該直線，故需要讓總的預測誤差

$J(beta_0, beta_1) = sum_{i=1}^n (y_i - hat{y}_i)^2 = sum_{i=1}^n [y_i - (beta_0+beta_1 x_i)]^2$

最小。

於是問題轉化為優化問題，選取 $beta_0, beta_1$ 使得

$min , J(beta_0, beta_1) = sum_{i=1}^n [y_i - (beta_0+beta_1 x_i)]^2$ (1)

這就是「最小二乘法」，有著很直觀的幾何解釋。

二、問題(1)求解

這是個求二元函數極小值問題。

根據微積分知識，二元函數極值是在一階偏導等於0點處取到：

$left{ begin{array}{l} dfrac{partial J}{partial beta_0} = -2 sumlimits_{i=1}^n big[ y_i - beta_0 - beta_1 x_i ) big] =0 [0.1cm] dfrac{partial J}{partial beta_1} = -2 sumlimits_{i=1}^n big[ y_i - beta_0 - beta_1 x_i ) big] x_i =0 end{array} right.$

解關於 $beta_0, beta_1$ 的二元一次方程組得

$left{ begin{array}{l} beta_0 = bar{y} - beta_1 bar{x} [0.2cm] beta_1 = dfrac{sumlimits_{i=1}^n x_i y_i - bar{y} sum_limits{i=1}^n x_i}{sumlimits_{i=1}^n x^2_i - bar{x} sumlimits_{i=1}^n x_i} = dfrac{sumlimits_{i=1}^n (x_i - bar{x})(y_i - bar{y})}{sumlimits_{i=1}^n (x_i - bar{x})^2} end{array} right.$

其中，

$left{ begin{array}{l} bar{x} = dfrac{1}{n} sumlimits_{i=1}^n x_i [0.2cm] bar{y} = dfrac{1}{n} sumlimits_{i=1}^n y_i end{array} right.$

三、提升：矩陣形式推導

將線性模型的全部預測值，用矩陣來寫：

$begin{bmatrix} 1 & x_1 vdots & vdots 1 & x_n end{bmatrix}_{n times 2} begin{bmatrix} beta_0 beta_1 end{bmatrix}_{2 times 1} =begin{bmatrix} hat{y}_1 vdots hat{y}_n end{bmatrix}_{n times 1}$

記

$X = begin{bmatrix} 1 & x_1 vdots & vdots 1 & x_n end{bmatrix}_{n times 2} , quad beta = begin{bmatrix} beta_0 beta_1 end{bmatrix}_{2 times 1} , quad hat{Y} =begin{bmatrix} hat{y}_1 vdots hat{y}_n end{bmatrix}_{n times 1}$

則矩陣表示為

$hat{Y}= X beta$

於是，讓預測誤差最小的「最小二乘法」優化問題就表示為

$minlimits_{beta} J(beta) = |Y - hat{Y} |^2 = |Y - X beta |^2$ (2)

這裡 $| cdot|$ 即向量的長度，計算一下：

$begin{eqnarray*} |Y - X beta |^2 &=& langle Y-X beta, Y-X beta rangle &=& (Y - X beta)^T (Y - X beta) &=& (Y^T - beta^T X^T)(Y - X beta) &=& Y^T Y - Y^T X beta - beta^T X^T Y + beta^T X^T X beta end{eqnarray*}$

同樣 $J(beta)$ 的極小值，在其一階偏導 $=0$ 處取到，為了對 $beta$ 求導，引入矩陣微積分知識：

$frac{partial x^T A}{partial x} = frac{partial A^T x}{partial x} = A$

$frac{partial x^T A x}{partial x} = A x + A^T x$

於是，令

$frac{partial |Y-X beta|^2}{partial beta}= frac{partial big( Y^T Y - Y^T X beta - beta^T X^T Y+beta^T X^T X beta big)}{partial beta} = 2X^T X beta - 2X^T Y = 0$

若 $X$ 滿秩，則 $X^T X$ 可逆，從而可得 $beta= (X^T X)^{-1} X^T Y$ .

注1：最小二乘法求解需要矩陣 $X$ 滿秩，否則只能採用梯度下降法求解.

注2：矩陣形式非常便於Matlab求解，同時也可以很容易地推廣到多元線性回歸：取

$X = begin{bmatrix} 1 & x_1^1 & cdots & x_m^1 vdots & vdots & ddots & vdots 1 & x_1^n & cdots & x_m^n end{bmatrix}$

即可，對應 $m$ 個自變數， $n$ 組樣本，第 $i$ 組樣本為 $(x_1^i, x_2^i, cdots, x_m^i, y_i)$ ，結果形式不變： $beta= (X^T X)^{-1} X^T Y$ , 此時 $beta = [beta_0, beta_1, cdots, beta_m]^T$ .

注3：一些非線性回歸也可以轉化為線性回歸來做：

例如，人口指數增長模型 $y=a e^{bx}$ ，做對數變換： $ln y = ln a + bx$ . 即將 $y$ 的數據取對數作為因變數，再與自變數 $x$ 數據做線性回歸，得到回歸係數 $beta_0, beta_1$ ，再由 $beta_0 = ln a, beta_1 = b$ 可得到 $a= e^{beta_0}, b= beta_1$ .

其它可變換為線性回歸的函數形式：

再例如，自變數 $x_1, x_2$ , 因變數 $y$ ，想做如下非線性回歸：

$y=beta_0 + beta_1 x_1 + beta_2 x_2 + beta_3 x_1^2 + beta_4 x_1 x_2 + beta_5 x_2^2$

令 $X_1 = x_1, X_2 = x_2, X_3 = x_1^2, X_4 = x_1 x_2, X_5= x_2^2$ ，做 $5$ 元線性回歸即可。

四、最小二乘法的概率解釋

真實數據中，一個 $x$ 值可能對應多個 $y$ 值，因為實際 $y$ 值可能受多種因素的影響，故可以假設任意一個 $x$ 值對應的 $y$ 的真實值是服從正態分布的隨機變數。

那麼，什麼時候擬合效果最好，當然是預測值 $hat{Y}=X beta$ 取值概率最大的時候，即最大似然原理。

記

$Y=Xbeta + varepsilon$

其中，

$Y = begin{bmatrix} y_1 vdots y_n end{bmatrix}$ 為真實值， $varepsilon = begin{bmatrix} varepsilon_1 vdots varepsilon_n end{bmatrix}$ 為預測誤差（即不能被線性模型刻畫的部分）。

假設 $varepsilon_i$ 為獨立同分布，服從均值為 $0$ ，方差為 $sigma^2$ 的正態分布。

該假設由中心極限定理可以保證，順便說一句，理想的誤差都得服從 $0$ 均值，等方差的正態分布，否則說明建立的模型不充分，誤差中還有趨勢沒有被提取出來。

則 $varepsilon_i$ 的概率密度為

$p(varepsilon_i)=frac{1}{sqrt{2 pi } sigma} exp Big( - frac{varepsilon_i^2}{2 sigma^2} Big), quad i = 1, cdots, n$

這就意味著，在給定 $x_i$ 和參數 $beta$ 情況下，因變數 $y_i$ 也服從正態分布，即

$p(y_i|x_i; beta)=frac{1}{sqrt{2 pi } sigma} exp Big( - frac{(y_i - x_i beta)^2}{2 sigma^2} Big), quad i=1, cdots, n$

其中， $x_i = (x_1^i, cdots, x_m^i)$

定義所有樣本數據關於參數 $beta$ 的似然函數：

$L(beta) = p(y|x;beta)=prod_{i=1}^n p(y_i | x_i ; beta) = prod_{i=1}^n frac{1}{sqrt{2 pi } sigma} exp Big( - frac{(y_i - x_i beta)^2}{2 sigma^2} Big)$

為了便於最大化似然函數，先取對數：

$ln big( L(beta) big) = ln bigg(prod_{i=1}^n frac{1}{sqrt{2 pi } sigma} exp Big( - frac{(y_i - x_i beta)^2}{2 sigma^2} Big) bigg)=n ln Big( frac{1}{sqrt{2 pi} sigma} Big) - sum_{i=1}^n frac{ (y_i - x_i beta)^2}{2 sigma^2}$