1. 程式人生 > >【數理統計基礎】 05 - 回歸分析

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

關於 以及 區間估計 否則 del 相互 不同之處 最小二乘 研究

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

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 - 回歸分析