更新時間：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\)稱為響應變量或因變量，\(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個樣本，實際上樣本也是一個隨機向量，那么對於第一層面，我們就可以改寫成：
將上述模型寫成矩陣的形式，其中\(X\)被稱為設計矩陣
2.2.1 Guass-Markov假設
因為我們知道隨機變量是可以求期望，求方差的。而在統計中，線性模型是基於諸多主觀假定的，其中我們常用的的假定是Guass-Markov假設。它主要假定了隨機誤差的期望是0，並且同方差，都是\(\sigma^2\)，此外隨機誤差之間都是無關的。
將Gauss-Markov假設寫成表達式的形式：
- \(D(\varepsilon_i) = \sigma^2\)
- \(cov(\varepsilon_i) = 0\)
- \(E(\varepsilon_i) = 0\)
其中，如果將假設1和2進行合並，我們就可以將假設簡化為
- \(cov(\varepsilon) = \sigma^2I_n\)
- \(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\)。因此，我們也把線性回歸稱為均值回歸.
其中，\(E(\varepsilon)\)為零向量
2.2.3 X固定下Y的分布
實際上，因為X是固定的，\(\beta\)是雖是未知參數，但也是固定的。因此, 有
說白了，就是\(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)\)
用圖來直觀的表示：
而當對\(\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 = X\beta + \varepsilon\)，可是可以通過實際的數值估計出\(\beta\)的值
2.4 其他
機器學習中有時也寫成向量的形式\(f(x) = w^Tx+b\)，將截距和其他回歸參數分開，但其實是沒有本質的區別的。
3. 基於損失函數的估計
下面主要是從隨機變量的層面出發來繼續深入討論線性模型，
在數學和計算機領域中，通常是使用損失函數來進行參數的估計，也可以說是模型的擬合。而值得注意的是，我們上述采用了的Guass-Markov假設，在參數\(\beta\)的估計是不需要用到的，這只是在我們之后的檢驗中需要使用到。但是我們不能說參數\(\beta\)的估計和我們給定的假設沒有一點關系，實際上根據不同的假設，我們選定的損失函數也是不同的。這一點，我會在第二篇“基於似然函數的估計”中提到。
先來談一下基於損失函數估計的原理：
- 定義模型中第\(i\)個殘差為：\(e_i = y_i - \hat y_i\)，整個殘差向量就可以寫成：\(e = Y - \hat Y\)，其中，\(\hat Y\)是擬合結果。
- 主觀選擇一個損失函數的形式：\(\rho(e_i)\)，又因為損失函數應該是包含待估參數\(\beta\)，因此損失函數又可以寫成是\(L(\beta)\)
- 利用損失函數最小化得到參數的估計：\(\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\)是一個數）：
要使損失函數最小化就是對損失函數求導，取到它的最小值（而這個損失函數是一個凸函數）
- 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函數的一個特例，它的圖像如下：
可以看見當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 無偏性
討論這個估計是否是無偏估計
4.3 有效性
經過上面的證明，我們可以得到：
其中，記\(H = X(X^TX)^{-1}X^T\)，且稱\(H\)為帽子矩陣，這個帽子矩陣有優良的性質，是一個冪等對稱矩陣。即有\(H^T = H\)、\(H^2 = H\)。此外，也把\(H\)稱為投影矩陣，即可以將隨機向量\(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\)的標准誤差為：
在這里，我們終於見到了假設中的\(\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\)，即有：
其中\(p\)為變量個數加上常數個數（即設計矩陣\(X\)的維數）
6. 估計量的區間估計
上面是做的估計都是點估計，下面談談怎么做區間估計
我們的目標是找到\(\bar \beta\)和\(\underline \beta\)，使得
為了實現這一目標，我們來構造含未知參數\(\beta\)的樞軸量
對於\(\beta_i\)，我們可以對它進行標准化，有：
此外，我們考慮\(\frac{\hat \sigma}{\sigma}\)
其中，\(\displaystyle \sum_{i=1}^n (\frac {y_i - \hat y_i}{\sigma})^2 \sim \chi^2(n-p)\)
將上述(I)和(II)結合起來，有
即有區間估計，\(\beta_i \underline + t_{1 - \frac{\alpha}{2}}(n-p) \hat{st.dev}(\hat \beta)\)
7. 結語
雖然最小二乘法的估計有明確的表達式，但是實際用計算機來擬合模型的時候，都不是直接用它估計的表達式來計算的，而是使用梯度下降的方法來加快收斂速度。在李航那本經典的小藍書，有提到到統計學習方法的三要素分別是：模型、策略和算法。對應這里實際上就是，模型是指線性模型，策略是指損失函數最小，而算法就是計算機實現這一模型的方法，也就是梯度下降算法。