深入理解線性模型(一)---基于損失函數的估計

更新時間:2019.10.31

1. 引言

??無論是統計學還是機器學習,我們最先接觸的模型(統計中的參數模型,機器學習中的有監督學習)都是線性模型,一個是因為它“簡單”,另一個是因為它是其他許多模型的一個衍生基礎,這也是為什么實際生活中雖然大多數都是非線性的,而我們還是要學習線性模型的原因。
我最初學習線性模型,覺得線性模型不就是一條直線嗎,這很簡單啊臉紅。然而,這畢竟是一堆天才研究了許多年才出來的東西,這肯定是沒有想象中的簡單,雖然它的思想確實簡單拍桌子。實際上,有很多模型都可以看做是線性的(即可以寫成\(Y = X\beta + \varepsilon\)的形式),雖然它們畫出來可能是一條曲線,例如像對數線性回歸\(log(y) = \beta_0 + \beta_1x\),多項式回歸\(y = \beta_0 + \beta_1x + \beta_2x^2\),還有回歸解釋變量中含有啞變量(定性變量)的情況,回歸樣條的情況等等都可以看做是線性模型,而這一類模型也被稱為廣義線性模型。
??對于線性模型我們通常采用最小二乘法來計算,當然這不是說這種方法是最好的。但基于歷史的原因,由于以前計算資源的匱乏,手算各種復雜的模型基本是行不通的,但是基于最小二乘法的線性模型的參數估計是有明確表達式的,并且具有優良的統計性質。
??最后,我將分成三篇以三種不同的角度(分別是損失函數、似然函數和貝葉斯估計)來分析線性模型的\(\beta\)\(\sigma\)這兩個參數的估計。而在本篇,是從損失函數的角度出發。

2. 從三個層面來看線性模型

??我們先從總體(理論),樣本(隨機變量)和數據三個層面來觀察一下一般的線性模型:

2.1 總體層面

??從理論層面上看(無數據)有:

\[Y = \beta_0 + \beta_1X_1 + \beta_2X_2 + \cdots + \beta_{p-1}X_{p-1} + \varepsilon \]

??其中,\(Y\)稱為響應變量或因變量,\(X_i\)為預測變量/解釋變量/自變量,\(\varepsilon\)為隨機誤差變量。在這里,我們可以把\((X_1, X_2, \cdots , X_{p-1}, Y)\)看做是一個總體,從理論層面上來觀測整一個模型。

2.2 樣本層面

??假設有\((X_{11}, X_{21}, \cdots , X_{p-1,1}, Y_1),\quad(X_{12}, X_{22}, \cdots , X_{p-1,2}, Y_2),\quad\cdots,\quad(X_{1n}, X_{2n}, \cdots , X_{p-1,n}, Y_n)\)等n個樣本,實際上樣本也是一個隨機向量,那么對于第一層面,我們就可以改寫成:

\[Y = \begin{pmatrix} Y_1\\ Y_2\\ \vdots\\ Y_n \end{pmatrix} ,\quad X = \begin{pmatrix} 1 & X_{11} & X_{21} & \cdots & X_{p-1,1}\\ 1 & X_{12} & X_{22} & \cdots & X_{p-1,2}\\ \vdots & \vdots & \vdots & & \vdots \\ 1 & X_{1n} & X_{2n} & \cdots & X_{p-1,n} \end{pmatrix} ,\quad \beta = \begin{pmatrix} \beta_0\\ \beta_1\\ \vdots\\ \beta_{p-1} \end{pmatrix} ,\quad \varepsilon = \begin{pmatrix} \varepsilon_1\\ \varepsilon_2\\ \vdots\\ \varepsilon_n\\ \end{pmatrix} \]

??將上述模型寫成矩陣的形式,其中\(X\)被稱為設計矩陣

\[Y_{n*1} = X_{n*p}\beta_{p*1} + \varepsilon_{n*1} \]

2.2.1 Guass-Markov假設

??因為我們知道隨機變量是可以求期望,求方差的。而在統計中,線性模型是基于諸多主觀假定的,其中我們常用的的假定是Guass-Markov假設。它主要假定了隨機誤差的期望是0,并且同方差,都是\(\sigma^2\),此外隨機誤差之間都是無關的。
??將Gauss-Markov假設寫成表達式的形式:

  1. \(D(\varepsilon_i) = \sigma^2\)
  2. \(cov(\varepsilon_i) = 0\)
  3. \(E(\varepsilon_i) = 0\)

??其中,如果將假設1和2進行合并,我們就可以將假設簡化為

  1. \(cov(\varepsilon) = \sigma^2I_n\)
  2. \(E(\varepsilon_i) = 0\)

??其中,\(cov(\varepsilon)\)是指隨機誤差向量的協方差陣,\(I_n\)\(n*n\)維的單位矩陣。此外,我們在做線性模型的時候,通常都是假設解釋變量\(X_i\)之間是無關的

2.2.2 均值向量

??當我們固定X(即相當于知道它的數值)的時候,結合Guass-Markov中\(E(\varepsilon_i) = 0\)的假設,我們可以求得固定\(X\)\(Y\)的條件期望應該為\(X\beta\),而實際上\(X\beta\)是我們的理論擬合值\(\hat Y\)。因此,我們也把線性回歸稱為均值回歸.

\[\begin{equation} \begin{cases} Y = X\beta + \varepsilon\\\\ \\\\ E(\varepsilon_i) = 0 \end{cases} => \quad E(Y|X) = E(X\beta + \varepsilon|X) = X\beta + E(\varepsilon) = X\beta \end{equation} \]

??其中,\(E(\varepsilon)\)為零向量

2.2.3 X固定下Y的分布

??實際上,因為X是固定的,\(\beta\)是雖是未知參數,但也是固定的。因此, 有

\[cov(Y) = cov(X\beta + \varepsilon) = cov(\varepsilon) = \sigma^2I_n \]

??說白了,就是\(Y\)滿足一個均值為\(X\beta\),方差為\(\sigma^2I_n\)的一個多維分布\(F_n(X\beta, \sigma^2I_n)\)\(Y_i\)滿足一個均值為\(X_i\beta\),方差為\(\sigma^2\)的分布\(F(X_i\beta, \sigma^2)\)

??用圖來直觀的表示:
固定x下y的分布圖

??而當對\(\varepsilon\)加上一個正態性的假設的時候,就能夠得到隨機誤差向量\(\varepsilon\)\(Y\)滿足多維正態分布\(\varepsilon \sim N(O, \sigma^2I_n)\)\(Y \sim N(X\beta, \sigma^2I_n)\)

2.3 數據層面

??當我們給定觀測值 \((x_{11}, x_{21}, \cdots , x_{p-1,1}, y_1),\quad(x_{12}, x_{22}, \cdots , x_{p-1,2}, y_2),\quad\cdots,\quad(x_{1n}, x_{2n}, \cdots , x_{p-1,n}, y_n)\),我們可以將\(Y\)\(X\)\(\beta\)\(\varepsilon\)表示成:

\[Y = \begin{pmatrix} y_1\\ y_2\\ \vdots\\ y_n \end{pmatrix} ,\quad X = \begin{pmatrix} 1 & x_{11} & x_{21} & \cdots & x_{p-1,1}\\ 1 & x_{12} & x_{22} & \cdots & x_{p-1,2}\\ \vdots & \vdots & \vdots & & \vdots \\ 1 & x_{1n} & x_{2n} & \cdots & x_{p-1,n} \end{pmatrix} ,\quad \beta = \begin{pmatrix} \beta_0\\ \beta_1\\ \vdots\\ \beta_{p-1} \end{pmatrix} ,\quad \varepsilon = \begin{pmatrix} \epsilon_1\\ \epsilon_2\\ \vdots\\ \epsilon_n\\ \end{pmatrix} \]

此時,模型仍可以寫成\(Y = X\beta + \varepsilon\),可是可以通過實際的數值估計出\(\beta\)的值

2.4 其他

??機器學習中有時也寫成向量的形式\(f(x) = w^Tx+b\),將截距和其他回歸參數分開,但其實是沒有本質的區別的。

3. 基于損失函數的估計

??下面主要是從隨機變量的層面出發來繼續深入討論線性模型,
??在數學和計算機領域中,通常是使用損失函數來進行參數的估計,也可以說是模型的擬合。而值得注意的是,我們上述采用了的Guass-Markov假設,在參數\(\beta\)的估計是不需要用到的,這只是在我們之后的檢驗中需要使用到。但是我們不能說參數\(\beta\)的估計和我們給定的假設沒有一點關系,實際上根據不同的假設,我們選定的損失函數也是不同的。這一點,我會在第二篇“基于似然函數的估計”中提到。
??先來談一下基于損失函數估計的原理:

  1. 定義模型中第\(i\)個殘差為:\(e_i = y_i - \hat y_i\),整個殘差向量就可以寫成:\(e = Y - \hat Y\),其中,\(\hat Y\)是擬合結果。
  2. 主觀選擇一個損失函數的形式:\(\rho(e_i)\),又因為損失函數應該是包含待估參數\(\beta\),因此損失函數又可以寫成是\(L(\beta)\)
  3. 利用損失函數最小化得到參數的估計:\(\hat \beta = \underset{\beta}{\arg \min} \hspace{1mm} \rho(e_i) = \underset {\beta}{\arg \min} \hspace{1mm} L(\beta)\),這一串長長的公式是說使得損失函數最小的變元\(\beta\)就是我們所要的估計參數\(\hat \beta\)

3.1 二次損失

??在做線性模型的時候,我們最常使用的是基于二次損失函數的最小二乘法,定義二次損失函數有:\(\sum_{i=1}^n e_i^2\) ,根據我們上面所說的均值回歸,以矩陣的形式有:
\begin{equation}
\begin{cases}
L(\beta) = \sum_{i=1}^n e_i^2 = \sum_{i=1}^n(y_i - \hat y_i)^2\\
\\
\hat Y = X\beta
\end{cases}
=> \sum_{i=1}^n (y_i - X_i\beta)^2 = (Y - X\beta)^T(Y - X\beta)
\end{equation}

??進一步化簡有(其中,\(Y^TX\beta\)是一個數):

\[(Y - X\beta)^T(Y - X\beta) = Y^TY - 2Y^TX\beta + \beta^TX^TX\beta \]

??要使損失函數最小化就是對損失函數求導,取到它的最小值(而這個損失函數是一個凸函數)

  • tips:值得注意的是,下文中有提到\(\hat Y = X\hat \beta\),嚴格來說上面提到的\(\hat Y\)是真實的擬合值,而下文提到的\(\hat Y\)是預測的擬合值

3.1.1 損失函數是最小的

??在求導之前,我們先思考一個問題,為什么\(\sum_{i=1}^n (y_i - X_i\beta)^2\)是最小的(即為什么不選擇其他\(\sum_{i=1}^n (y_i - ?)^2\),為什么前面的要小于等于\(\sum_{i=1}^n (y_i - ?)^2\))。實際上,這就相當于對于\(E(X - ?)^2\)來說,當?取什么時,這個式子達到最小。而我們知道,當?\(X\)的期望\(E(X)\),這個式子就能達到最小。同理(求和和期望只是差了一個系數),對于\(\sum_{i=1}^n (y_i - ?)^2\))來說,當?\(Y_i\)的期望\(E(Y_i)\),這個式子就能達到最小,而這個\(E(Y_i)\)恰好是\(X_i\beta\),也就是\(\hat Y_i\)實際上是這組數據中,\(Y\)的均值。這樣也就保證我們選擇的損失函數是在這種二次形式中是最好。

3.1.2 損失函數最小化

??下面對損失函數進行求導:
\begin{equation}
\begin{split}
\frac {\mathrm{d}L(\beta)}{\mathrm{d}\beta} & = \frac {\mathrm{d} Y^TY - 2Y^TX\beta + \beta^TX\beta}{\mathrm{d} \beta}\\
& = 0 - 2X^TY + 2X^TX\beta
\end{split}
\end{equation}

??令\(\frac {\mathrm{d}L(\beta)}{\mathrm{d}\beta} = 0\),有:

\begin{equation}
0 - 2X^T Y + 2X^TX\beta = 0 => \hat \beta =
\begin{cases}
(X^T X)^{-1} X^TY, & (X^TX)可逆 \\
\\
(X^T X)^- X^TY, & (X^TX)不可逆
\end{cases}
\end{equation}

??正如之前提到的,我們通常假定解釋變量\(X_i\)是無關的(即列滿秩),此時便有\((X^TX)\)可逆(因為\(Ax=0\)\(A^TAx=0\)是同解方程組,所以可以證出\(r(A)=r(A^TA)\)的),所以通常將待估參數表示為\(\hat \beta = (X^TX)^{-1}X^TY\)
??從上述的步驟中,我們是可以看到做估計的時候沒有用到Guass-Markov的假設。而在數學和機器學習中,有時也就直接使用參數估計后的模型進行測試,然后求出它的均方誤差,并沒有對這個估計量做檢驗。而實際上,\(\hat \beta\)最佳的線性無偏估計,說它最佳是因為它方差是最小的(檢驗估計量的有效性中將會提到)。
??最后值得一提的是,上面所述的二次損失形式其實是基于方差齊性\(cov(\varepsilon)=\sigma^2I_n\)的假設,當我們改變這個假設的時候,二次損失的將會有所不同,這個我會在第二篇“基于似然函數的估計”中會進一步提到。

3.2 其他損失函數

??既然提到了二次損失,在這里也提一下其他的一些損失函數

3.2.1 最小絕對損失

??定義為:\(L(\beta) = \sum_{i=1}^n |y_i - x_i\beta|\),也叫一次損失
??此外,對于\(E(|y_i - ?|)\),當?取中位數的時候,這個損失函數是最小的,所以有時也稱為最小中位估計。而由于中位數比平均數要穩健,所以對比二次損失,絕對損失對異常點不敏感。當然這里,當我們假定\(\varepsilon\)是服從均值為0的正態分布的時候,\(X\beta\)實際上也是固定X下Y的中位數

3.2.2 huber函數

??huber函數主要是用于穩健回歸的,其中絕對損失其實就是huber函數的一個特例,它的圖像如下:
huber函數圖
??可以看見當u到一定程度的時候,它的損失是呈水平的,而異常點是指那些遠離理論擬合值的觀測點,即指u很大的點,這一說明了huber函數的穩健(對異常點不敏感)。當然huber函數也存在一定的缺陷,像多一個參數k,還有不光滑的部分。

3.2.3 分位回歸的損失

??定義為:\(L(\beta) = \sum_{i=1}^n e_i(\tau - I(e_i < 0))\),其中I()是指示函數,\(0 < \tau < 1\)。上面所述的都是對稱的損失函數,但是分位回歸的損失函數是不對稱的,實際上不對稱帶來的影響會左右曲線的擬合,這時就需要使用分位回歸的損失。
??正如下面這個圖所示:
不對稱帶來的影響
??實際上,絕對損失同樣也可以看做是分位回歸的損失函數的一個特例,當\(\tau\)取0.5時,兩者實際上就差了一個常數。

4. 估計量\(\hat \beta\)的檢驗

??談過估計之后,繼續談一下估計量\(\hat \beta\)的檢驗,而實際上這個估計是最佳線性無偏估計

4.1 估計量的評價標準

??估計量的評價標準主要有三個,無偏性、有效性、一致性。以這里的\(\hat \beta\)為例,

  • 無偏性指的是估計量的期望等于參數的真實值,即\(E(\hat \beta) = \beta\),其中\(\beta\)實際上是一個隨機向量,因為X是固定的,而Y是隨機向量。
  • 有效性是指在諸多無偏估計中,方差最小的那一個估計量。
  • 一致性是描述當樣本容量無窮時,估計量限接近參數的真實值,即\(n-> + \infty, \hat \beta -> \beta\),表示成數學語言就是\(P \underset{n-> + \infty}(|\hat \beta - \beta| < \varepsilon) = 1\)
    ??在這之中,一致性是最容易滿足的,其次是無偏性,最后才是有效性。在這里,只討論無偏性和有效性

4.2 無偏性

??討論這個估計是否是無偏估計

\[E(\hat \beta) = E((X^TX)^{-1}X^TY) = (X^TX)^{-1}X^TE(Y) = (X^TX)^{-1}X^TX\beta = \beta \]

4.3 有效性

??經過上面的證明,我們可以得到:

\[ \hat Y = X \hat \beta = X(X^TX)^{-1}X^TY = HY \]

??其中,記\(H = X(X^TX)^{-1}X^T\),且稱\(H\)為帽子矩陣,這個帽子矩陣有優良的性質,是一個冪等對稱矩陣。即有\(H^T = H\)\(H^2 = H\)。此外,也把\(H\)稱為投影矩陣,即可以將隨機向量\(Y\)投影到\(X\)的超平面上。
??再來看一幅圖
Y在X平面的投影
??我們知道當垂直時,線段到平面的距離是最短的。即當\(HY\perp(I_n - H)Y\)時,殘差\(e\)是最小的,這也直觀地證明這個估計是方差最小的(因為殘差\(e\)是方差的一個度量),因此我們也稱這個估計\(\hat \beta\)一切線性無偏估計中方差最小的。
??既然\(\hat \beta\)是最小的,那么下面來推導下它的方差究竟是多少
\begin{equation}
\begin{split}
cov(\hat \beta) & = cov((X^T X)^{-1} X^T Y)\\
& = (X^T X)^{-1} X^T cov(Y) X (X^T X)^{-1}\\
& = (X^T X)^{-1} X^T \sigma^2 I_n X (X^T X)^{-1}\\
& = \sigma^2 (X^T X)^{-1}
\end{split}
\end{equation}

??因此,根據\(\hat \beta_i\)的協方差,我們可以得到\(\hat \beta_i\)的標準誤差為:

\[st.dev(\hat \beta_i) = \sigma \sqrt{(X^TX)^{-1}_{ii}} \]

??在這里,我們終于見到了假設中的\(\sigma\)(同方差),而這正也說明了我之前所說的估計\(\beta\)不用假設,但是檢驗的時候就需要用到

5. \(\sigma^2\)的估計

??雖然算出了\(\hat \beta_i\)的標準誤差,但是它表達式中的\(\sigma\)實際上是一個未知參數,下面對其進行估計。我們知道方差是指偏離平均水平的程度,因此我們很自然地想到使用\(\frac{\displaystyle \sum_{i=1}^n(y_i - \hat y_i)^2}{n}\)來估計,即使用\(\frac{\displaystyle \sum_{i=1}^n(e_i)^2}{n}\)來估計,因為我們之前說過\(\hat Y\)是X固定下Y的平均水平。因此,有
\begin{equation}
\begin{split}
\displaystyle \sum_{i=1}^n (e_i)^2 & = e^T e\\
& = Y^T(I_n - H)^T(I_n - H)Y \\
& = Y^T(I_n - H)Y
\end{split}
\end{equation}

??其中,由于\(H\)是冪等對稱陣,所以\((I_n - H)\)實際上也是冪等對稱陣。
??下面給出兩種證法:
??證法1:
\begin{equation}
\begin{split}
E(Y^T(I_n - H)Y) & = E(tr(Y^T (I_n - H)Y))\\
& = E[tr(YY^T (I_n - H))]\\
& = tr[(I_n - H)E(YY^T)]\\
& = tr(I_n - H)[cov(Y) + E(Y)E(Y)^T]\\
& = tr(I_n - H)[\sigma^2I_n + X\beta \beta^T XT]\\
& = (n-p)\sigma^2 + tr[(I_n - H)X\beta \beta^T XT]\\
& = (n-p)\sigma^2 + tr[(I_n - X(X^T X)^{-1} X^T)X\beta \beta^T XT]\\
& = (n-p)\sigma^2
\end{split}
\end{equation}

??證法2:
??\begin{equation}
\begin{split}
E(Y^T(I_n - H)Y) & = E(tr(Y^T(I_n - H)Y))\\
& = E(Y^T)(I_n - H)E(Y) + tr((I_n - H)cov(Y))\\
& = \beta^T X^T(I_n - H)X\beta + (n-p)\sigma^2\\
& = \beta^TX ^T(I_n - X(X^T X)^{-1} X^T)X\beta + (n-p)\sigma^2\\
& = (n-p)\sigma^2
\end{split}
\end{equation}

??其中,由于帽子矩陣的是冪等陣,所以特征值只有0或者1,設有p個特征值為1,因此\(I_n - H\)有n-p個特征值為1,即\(tr(I_n - H) = n-p\)。實際上,當X列滿秩的時候,有\(tr(H) = tr(X(X^TX)^{-1}X^T) = tr(X^TX(X^TX)^{-1}) = tr(I_p) = p\),因此可以看出H的秩等于\(X\)的維數,也等于\(X\)的變量個數加上常數。在這里,我們也更加清楚,為什么我們要假設解釋變量\(X_i\)都是無關的。
??由上面的推導過程可以看出,\(\sigma^2\)的無偏估計應該是\(\frac {\displaystyle \sum_{i=1}^n e_i}{n-p}\),將殘差平方和\(\displaystyle \sum_{i=1}^n e_i\)記為\(SSE\),并且記均方誤差為\(MSE\),即有:

\[\hat \sigma^2 = \frac {SSE}{n-p} = MSE \]

??其中\(p\)為變量個數加上常數個數(即設計矩陣\(X\)的維數)

6. 估計量的區間估計

??上面是做的估計都是點估計,下面談談怎么做區間估計
??我們的目標是找到\(\bar \beta\)\(\underline \beta\),使得

\[P(\bar \beta_i < \beta_i < \underline \beta_i) = 1 - \alpha \]

??為了實現這一目標,我們來構造含未知參數\(\beta\)的樞軸量
??對于\(\beta_i\),我們可以對它進行標準化,有:

\[\frac {\hat \beta_i - \beta_i}{\sigma \sqrt{(X^TX)^{-1}_{ii}}} \sim N(0, 1) \tag{I} \]

??此外,我們考慮\(\frac{\hat \sigma}{\sigma}\)

\[ \frac{\hat \sigma}{\sigma} = \sqrt{\frac {\frac {\displaystyle \sum_{i=1}^n (y_i - \hat y_i)^2}{n-p}}{\sigma^2}} = \sqrt{\frac {\displaystyle \sum_{i=1}^n (\frac {y_i - \hat y_i}{\sigma})^2}{n-p}} \tag{II} \]

??其中,\(\displaystyle \sum_{i=1}^n (\frac {y_i - \hat y_i}{\sigma})^2 \sim \chi^2(n-p)\)
??將上述(I)和(II)結合起來,有

\[\frac {\hat \beta_i - \beta_i}{\hat \sigma \sqrt{(X^TX)^{-1}_{ii}}} = \frac{\frac {\hat \beta_i - \beta_i}{\sigma \sqrt{(X^TX)^{-1}_{ii}}}} {\frac {\hat \sigma}{\sigma}} = \frac{\frac {\hat \beta_i - \beta_i}{\sigma \sqrt{(X^TX)^{-1}_{ii}}}} {\sqrt{\frac {\displaystyle \sum_{i=1}^n (\frac {y_i - \hat y_i}{\sigma})^2}{n-p}}} \sim t(n-p) \]

??即有區間估計,\(\beta_i \underline + t_{1 - \frac{\alpha}{2}}(n-p) \hat{st.dev}(\hat \beta)\)

7. 結語

??雖然最小二乘法的估計有明確的表達式,但是實際用計算機來擬合模型的時候,都不是直接用它估計的表達式來計算的,而是使用梯度下降的方法來加快收斂速度。在李航那本經典的小藍書,有提到到統計學習方法的三要素分別是:模型、策略和算法。對應這里實際上就是,模型是指線性模型,策略是指損失函數最小,而算法就是計算機實現這一模型的方法,也就是梯度下降算法。

posted @ 2019-10-31 16:19  jianli-Alex  閱讀(...)  評論(...編輯  收藏
11选5走势图