參數估計和假設檢驗是數理統計的兩個基礎問題，它們不光運用於常見的分布，還會出現在各種問題的討論中。本篇開始研究另一大類問題，就是討論多個隨機變量之間的關系。現實生活中的數據雜亂無章，夠挖掘出各種變量之間的關系非常有用，它可以預估變量的走勢，能幫助分析狀態的根源。關系分析的著手點可以有很多，我們從最簡單直觀的開始，逐步展開討論。

1. 一元線性回歸

1.1 回歸分析

　　如果把每個量都當做隨機變量，問題的討論會比較困難，或者得到的結論會比較受限。一個明智做法就是只把待考察的量\(Y\)看做隨機變量，而把其它量\(X_i\)看成是自主選定的。即使都看成變量，也是把\(Y\)看成因變量，而把\(X_i\)看成自變量。該模型同樣是研究某個隨機變量的情況，不同之處在於更加關註變量與各因素的函數關系，希望能找到影響隨機變量的主要因素並給出表達式。

　　如式（1）所示，選定要關註的因素\(X_i\)，並假定它們以函數\(f(X_1,\cdots,X_p)\)形式影響變量\(Y\)，其它的因素統一放到隨機變量\(e\)中。其中函數\(f\)稱為\(Y\)對\(X_i\)的回歸函數或回歸方程，\(e\)則是隨機誤差。由於已經提取出主要因素，這裏假定\(e\)的均值為\(0\)，並且它是與\(f\)獨立的。在應用場景，一般給定回歸函數一個含參表達式（比如後面的線性回歸），這樣的問題稱為參數回歸問題，否則叫非參數回歸問題。

\[Y=f(X_1,X_2,\cdots,X_p)+e,\;\;E(e)=0,\;0<D(e)<\infty\tag{1}\]

　　在微積分中我們知道，多元函數\(f(x_1,\cdots,x_p)\)一般可以在任何點進行泰勒展開，其中最簡單的就是線性展開。線性關系由於其形式簡單，以及在局部能很好地逼近函數，在數學的各分支都被重點討論。在回歸分析中，這樣的模型便稱為線性回歸，這裏先從最簡單的一元線性回歸討論起。

　　一元線性回歸的模型是式（2）左，在提取出線性關系\(b_0+b_1X\)後，\(Y\)的剩余因素或隨機性就都落在隨機變量\(e\)上。所以從另一個角度看，回歸分析是要找出隨機變量的“確定”部分和隨機部分，這種分解更能幫助分析隨機現象。自然地，分析是基於\(n\)個樣本點\((X_i,Y_i)\)，其中\(X_i\)可能也是隨機產生的，但在這個模型裏一律看做定量。還要註意，這時每個\(Y_i\)是\(e_i\)的一個位移，它不再與\(Y\)同分布。

\[Y=b_0+b_1X+e;\;\;Y_i=b_0+b_1X_i+e_i=a_0+a_1(X_i-\bar{X})+e_i\tag{2}\]

1.2 系數點估計

　　每一次試驗相互獨立，因此得到\(n\)個獨立變量（式（2）右），其中\(b_0,b_1\)是待定系數。為了方便計算和討論，一般還會把式（2）右的線性部分“中心化”。問題等價於討論\(a_0,a_1\)的值，但要註意這裏\(\bar{X}\)是依賴於具體樣本的。在得知樣本點\((X_i,Y_i)\)的情況下，如何確定系數比較合理？在式（2）中，我們把\(Y_i\)看做是有誤差\(e_i\)的變量，因此讓誤差的平方和達到最小是一個比較好的模型。

　　式（3）取最小值時\(a_0,a_1\)便是合理的參數估計，利用偏導為零容易算得式（4）中對\(a_0,a_1\)的估計，這個結論非常得益於剛才的中心化。求解的方法其實就是最小二乘法，這在後面再展開討論。\(a_0\)表示\(Y\)的中心，估計值\(\alpha_0\)十分合理。\(a_1\)應當是\(Y\)關於\(X\)的斜率，單點的斜率是\(\dfrac{Y_i-\bar{Y}}{X_i-\bar{X}}\)，將分子分母同時乘以\(X_i-\bar{X}\)並相加，化簡後便得到\(\alpha_1\)，故它也是斜率的合理估計。

\[\sum\limits_{i=1}^n(Y_i-a_0-b_1X_i)^2\tag{3}\]

\[\alpha_0=\bar{Y};\;\;\alpha_1=\sum_{i=1}^n\dfrac{X_i-\bar{X}}{S^2}Y_i,\;S^2=\sum_{i=1}^n(X_i-\bar{X})^2\tag{4}\]

　　另外還要註意，式（4）中\(X_i\)是定值，而\(Y_i\)獨立隨機變量，\(\alpha_0,\alpha_1\)都是\(Y_i\)的線性函數，這對於下面的討論很重要。估計合理的另一個基本要求應當是誤差估計、即統計量（隨機變量）\(\alpha_0,\alpha_1\)的期望值應當就是\(a_0,a_1\)，利用式（5）左很容易驗證結論成立（式（5）右）。以下令\(e\)的方差為\(\sigma^2\)，利用\(Y_i\)的無關性也容易有式（6）。其中\(D(\alpha_1)\)分母中的\(S^2\)有直觀的含義，當\(X_i\)比較分散時，得到的斜率估計越準確。另外還可以證明，\(\alpha_0,\alpha_1\)是\(a_0,a_1\)的MVU估計。

\[E(Y_i)=a_0+a_1(X_i-\bar{X})\;\Rightarrow\;E(\alpha_0)=a_0,\;E(\alpha_1)=a_1\tag{5}\]

\[D(Y_i)=\sigma^2\;\Rightarrow\;D(\alpha_0)=\dfrac{\sigma^2}{n},\;D(\alpha_1)=\dfrac{\sigma^2}{S^2}\tag{6}\]

　　還有一點，把\(\alpha_0,\alpha_1\)看成\(Y_i\)的線性函數，觀察兩者的“系數向量”，發現它們內積為\(0\)。從向量的角度它們就是直交的，經驗證\(\alpha_0,\alpha_1\)也的確是（線性）不相關的，這個結論非常重要，也顯示了前面中心化的意義。另外，當\(e\)是正態分布時，\(\alpha_0,\alpha_1\)也都是正態分布，故可知它們獨立。

1.3 誤差估計

　　對於模型（2）來說，目前還有\(e\)的方差\(\sigma^2\)沒有討論，在有了系數估計（4）後，現在來估計誤差的方差。隨著\(X\)的變化，\(Y\)的中心也跟著變化，其誤差的方差自然也要以具體的中心為準。在樣本點\((X_i,Y_i)\)處，誤差\(\delta_i\)（式（7））也稱為殘差，它們的平均平方和理應作為方差的估計。但由於\(\alpha_0,\alpha_1\)的估計中消耗了兩個自由度，故可驗證式（7）才是\(\sigma^2\)的無偏估計。

\[\hat{\sigma}^2=\dfrac{1}{n-2}\sum_{i=1}^n\delta_i^2,\;\;\delta_i=Y_i-\alpha_0-\alpha_1(X_i-\bar{X})\tag{7}\]

　　具體計算步驟參考教材（或自行證明），結果是得到式（8），這樣就不難得到是（7）了。當然在實際計算時，可以直接展開得到式（9），然後利用現成的\(X_i,Y_i,\alpha_j\)來加速計算。而且從式（9）中還能得到更有用的結論，註意其中的後兩項\(n\bar{Y}^2=Z_1^2\)和\(S^2\alpha_1^2=Z_2^2\)，\(Z_1,Z_2\)都是\(Y_i\)的線性函數，且系數向量是兩個相互正交的標準化向量。

\[\sum_{i=1}^n\delta_i^2=\sum_{i=1}^n(e_i-\bar{e})^2+\dfrac{1}{S^2}\left(\sum_{i=1}^n(X_i-\bar{X})e_i\right)^2\tag{8}\]

\[\sum_{i=1}^n\delta_i^2=\sum_{i=1}^nY_i^2-n\bar{Y}^2-S^2\alpha_1^2\tag{9}\]

　　當\(e\)是正態分布時，\(Y_i\)也是正態分布，利用正交變換的性質，易知式（9）等於\(Z_3^2+\cdots+Z_n^2\)，其中\(Z_j\sim N(0,\sigma^2)\)，這便容易有式（10）的結論。關於殘差，還有兩點需要註意，式（8）如果很大或者殘差體現出某些規律性，則說明線性模型不太合適，或還有重要因素沒有被提取出來。

\[e\sim N(0,\sigma^2)\;\Rightarrow\;\sum_{i=1}^n\dfrac{\delta_i^2}{\sigma^2}\sim\chi_{n-2}^2\tag{10}\]

1.4 區間估計

　　有了點估計，便可以做區間估計，為了能使用樞軸函數，這裏還是假定\(e\)為正態分布。首先由公式（5）（6）可知\(\alpha_0,\alpha_1\)滿足式（11）的分布，當\(\sigma\)已知時，樞軸函數很容易得到。當\(\sigma\)未知時，由剛才的討論知\(\alpha_0,\alpha_1\)與\(\hat{\sigma}^2\)是相互獨立的，這樣便能用\(\hat{\sigma}\)替代\(\sigma\)，得到式（12）的樞軸變量。

\[\alpha_0\sim N(a_0,\dfrac{\sigma^2}{n});\;\;\alpha_1\sim N(a_1,\dfrac{\sigma^2}{S^2})\tag{11}\]

\[\dfrac{\sqrt{n}(\alpha_0-a_0)}{\hat{\sigma}}\sim t_{n-2};\;\;\dfrac{S(\alpha_1-a_1)}{\hat{\sigma}}\sim t_{n-2}\tag{12}\]

　　線性回歸的目的自然是為了進行預測，但在僅知道樣本點且把\(X_i\)看成定量的情況下，其實是無法估計最初式（2）左中的\(b_0,b_1\)的。因此要註意，在用\(y=a_0+a_1(x-\bar{X})\)預估\(Y\)時，我們不光丟失了誤差\(e\)，還丟失了\(X\)非連續得來的誤差。前者通過合理建模來降低誤差，後者則只能通過增加\(X_i\)的數量和密度來降低誤差。

　　這一點容易通過估計值\(y\)的方差看出（式（13））。首先在樣本數不變的情況下，\(x\)離樣本中心\(\bar{X}\)越近方差越小，這個結論符合直覺，樣本離預測點越近精度越高。另一方面，樣本數越大方差也越小，這個很好理解。結合這兩方面，當\(n\)足夠大且\(x\)離樣本中心足夠近，估計的方差就可以任意小。

\[D(y)=\left(\dfrac{1}{n}+\dfrac{(x-\bar{X})^2}{S^2}\right)\sigma^2\tag{13}\]

2. 多元線性回歸

2.1 系數估計

　　現實中的因變量可能受多個因素的影響，這些因素可能有主次之分，也可能是聯合作用。無論如何，為了對因變量進行更加深入細致的分析，必須加入更多的自變量進行分析。另外同樣的道理，多元函數在局部都可以用線性函數很好地近似，因此我們也可以建立式（14）中的模型和中心化樣本表達式。為表達方便，本段下面就直接把\(X_{ki}-\bar{X}_k\)記作\(X_{ki}\)。

\[Y=b_0+\sum_{k=1}^p b_kX_k+e;\;\;Y_i=a_0+\sum_{k=1}^p a_k(X_{ki}-\bar{X}_k)+e_i\tag{14}\]

　　多元模型的討論內容和方法與一元的差別不大，但直接的討論會很繁瑣，必須借助於線性代數的工具，請註意前後對比。為討論方便，首先規定式（15）的簡寫，並記\(\gamma\)的點估計為\(\alpha\)。然後定義列向量的期望\(E(\alpha)=[E(\alpha_i)]\)，以及協方差\(\text{Cov}(\alpha,\beta)=[\text{Cov}(\alpha_i,\beta_j)]\)，且不難驗證有式（16）成立。其實利用算子理論證明會很簡單，但光憑形式化的假設，也不難完成證明，請獨立嘗試。

\[\beta=\begin{bmatrix}Y_1\\Y_2\\\vdots\\Y_n\end{bmatrix},\;\gamma=\begin{bmatrix}a_0\\a_1\\\vdots\\a_p\end{bmatrix},\;A=\begin{bmatrix}1&\cdots&1\\X_{11}&\cdots&X_{1n}\\\vdots&\ddots&\vdots\\X_{p1}&\cdots&X_{pn}\end{bmatrix}\tag{15}\]

\[E(A\alpha)=AE(\alpha);\;\;\text{Cov}(A\alpha,B\beta)=A\text{Cov}(\alpha,\beta)B^T\tag{16}\]

　　有了矩陣的定義，就可以直接利用線性代數中最小二乘的結論，得到（17）左的正則方程，以及式（17）的\(\gamma\)最小二乘解。式（18）推導出\(\alpha_i\)是\(a_i\)的無偏估計，且都是\(Y_i\)的線性函數，繼而還可以得到式（19）的協方差公式。註意到\(A\)中除第一列外，每行的和都是\(0\)，故\(L\)及\(L^{-1}\)都具有形式\(\begin{bmatrix}1&0\\0&B_{p\times p}\end{bmatrix}\)。這說明\(\alpha_0=\bar{Y}\)與其它\(\alpha_i\)互不相關，這與一元的情況是一致的。

\[L\alpha=A\beta\;\Rightarrow\;\alpha=L^{-1}A\beta,\;\;(L=AA^T)\tag{17}\]

\[E(\alpha)=L^{-1}AE(\beta)=L^{-1}AA^T\gamma=\gamma\tag{18}\]

\[\text{Cov}(\alpha,\alpha)=L^{-1}A\text{Cov}(\beta,\beta)A^TL^{-1}=\sigma^2L^{-1}\tag{19}\]

2.2 誤差估計

　　現在來分析誤差\(e\)，首先記殘差向量\(\delta=\beta-A^T\alpha\)，容易證明\(E(\delta)=0\)，而且根據式（20）的推導可知\(\alpha_i\)與\(\delta_j\)互不相關。另外可以算得式（21）的協方差，其中\(B\)是一個秩為\(p+\)的非負定方陣。根據\(B^2=B\)可以證明，\(B\)的\(p+1\)個特征值都是\(1\)，從而它的跡\(tr(B)=p+1\)（主對角線之和，請參考線性代數）。

\[\text{Cov}(\hat{\alpha},\delta)=L^{-1}A\text{Cov}(\beta,\beta)(I_n-A^TL^{-1}A)=0\tag{20}\]

\[\text{Cov}(\delta,\delta)=(I_n-B)\text{Cov}(\beta,\beta)(I_n-B)=\sigma^2(I_n-B),\;\;(B=A^TL^{-1}A)\tag{21}\]

　　為了估計\(\sigma^2\)，自然想到討論殘差平方和\(\sum\limits_{i=1}^n\delta_i^2\)。式（22）計算了它的期望值，這樣就可以用式（23）來無偏估計\(\sigma^2\)。

\[E(\sum_{i=1}^n\delta_i^2)=\sum_{i=1}^nD(\delta_i)=tr(\text{Cov}(\delta,\delta))=\sigma^2(n-p-1)\tag{22}\]

\[\hat{\sigma}^2=\dfrac{1}{n-p-1}\sum_{i=1}^n\delta_i^2\tag{23}\]

　　殘差平方和\(\delta^T\delta\)是一個半正定二次型，展開整理後有式（24）成立，它滿足柯赫倫定理的條件。故假定\(e\)為正態分布的情況下，有式（25）左成立。另外由於\(\alpha_i\)與\(\delta_j\)互不相關，則\(\alpha_i\)與\(\hat{\sigma}\)也不相關，正態分布下它們還是相互獨立的，這就得到式（25）右的樞軸變量。

\[\beta^T\beta=\beta^TB\beta+\delta^T\delta\tag{24}\]

\[\sum_{i=1}^n\dfrac{\delta_i^2}{\sigma^2}\sim\chi_{n-p-1};\;\;\dfrac{\alpha_i-a_i}{\sqrt{L_{ii}^{-1}}\hat{\sigma}}\sim t_{n-p+1}\tag{25}\]

2.3 假設檢驗

　　線性回歸的假設往往是針對線性系數\(a_k\)的，如果僅是對單個系數的檢驗，直接利用式（25）的樞軸變量即可。實際應用中最常用的假設是\(a_k=0\)，它說明因素\(X_k\)對\(Y\)其實是不相關的，這對檢驗變量相關性很有用（但更偏重\(X_k\)對\(Y\)的影響）。觀察式（17），你會發現\(\alpha_k\)並不只與\(X_k,Y\)有關，它與上面的一次模型得到的結論不一樣。可以這樣解釋：更多因素的加入使得模型更加精確。

　　但是不是因素越多越好？如果加入的是真正影響\(Y\)的元素，對模型自然是有益的，否則多加入的元素只能增加隨機性，從而對結論精度造成影響。樣本不足的情況下，以上模型容易把無效元素估計成“假”的關系，從而影響真實因素的作用。但逐個地檢驗無效元素，有時效果並不好，因為元素之間的復雜關系和隨機性會使得檢驗出現較大偏差。

　　檢驗較多無關參數時，最好能將它們捆綁操作，當選定好要檢驗的無關參數後，甚至可以將將模型中的其它參數去除，以簡化討論，也就是說假設條件變成\(a_1=\cdots=a_p=0\)。但這個多變量的假設很難建立之前單變量的樞軸變量，我們需要另外找一個變量作為度量的對象。在鑒別“有效、無效”元素的問題中，註意“有效”的元素的典型特征，就是使得殘差平方和變小，或者說使得\(\hat{\sigma}^2\)盡量小。這便是我們要找的“值”，具體來說，就是要度量\(\hat{\sigma}^2\)和\(\sigma^2\)之間的差別。

　　但由於\(\sigma^2\)未知，必須找統計量來替代它，在假設條件下，自然是用\(S_Y^2\)。當直接用\(\hat{\sigma}^2\)和\(S_Y^2\)難產生好的樞軸變量，原因主要是系數的影響，這時我們自然想到直接比較殘差平方和。為此記\(R_1=\delta^T\delta\)，並記假設條件下的殘差平方和為\(R_2\)。為了討論方便，這裏把式（14）稍作修改，就是先作出估計\(\alpha_0=\bar{Y}\)，然後用\(Y_i\)取代\(Y_i-\alpha_0\)重新建模，隨之\(\gamma,A\)中的第一列也去除。

　　但新的模型仍然能得到式（17）的估計式，以及殘差向量\(\delta=\beta-A^T\alpha\)。這個模型下\(R_1,R_2\)如式（26）所示，不難發現\(R_2-R_1=\beta^TB\beta\)，而已知\(B^2=B\)，所以\(R_2-R_1\)也是一個半正定二次型。再次使用柯赫倫定理可有式（27）左成立，並且\(R_2-R_1\)與\(R_1\)互相獨立，這等價於與\(\hat{\sigma}^2\)互相獨立，所以得到式（27）右的樞軸變量。註意到\(R_2\geqslant R_1\)，故檢驗否定的條件應當是\(F>C\)。

\[R_1=\delta^T\delta=\beta^T(I_n-B)\beta;\;\;R_2=\beta_T\beta\tag{26}\]

\[\dfrac{R_2-R_1}{\sigma^2}\sim\chi_p^2;\;\;\dfrac{R_2-R_1}{r\hat{\sigma}^2}\sim F_{p,n-p-1}\tag{27}\]

【數理統計基礎】 05 - 回歸分析