數理統計2：為什么是正態分布，正態分布均值與方差的估計，卡方分布

上一篇文章提到了一大堆的統計量，但是沒有說到它們的用處。今天，我們就會接觸到部分估計量，進入到數理統計的第一大范疇——參數估計，同時也會開始使用R語言進行模擬。由於本系列為我獨自完成的，缺少審閱，如果有任何錯誤，歡迎在評論區中指出，謝謝！

目錄

• Part 1：為什么是正態分布
• Part 2：正態分布均值估計
• Part 3：正態分布方差估計
• Part 4：卡方分布

Part 1：為什么是正態分布

為什么要突然提到正態分布的參數估計？原因有以下幾個。首先，正態分布是生活中最常見的分布，許多隨機事件的分布可以用正態分布來概括。林德貝格勒維中心極限定理告訴我們，二階矩存在的獨立同分布隨機變量列\(\{\xi_n\}\)，記它們的和為\(S_n\)，\(\mathbb{E}(\xi_1)=\mu\)，\(\mathbb{D}(\xi_n)=\sigma^2\)，則

\[\frac{S_n-n\mu}{\sqrt{n}\sigma}\stackrel{d}\to N(0,1). \]

剛剛學完概率論的同學應該對這個結論不陌生。

而中心極限定理的條件實際上並不需要這么強，林德貝格費勒定理去除了同分布的約束，只要\(\{\xi_n\}\)滿足\(\forall \tau>0\)，

\[\frac{1}{\sum_{k=1}^n\mathbb{D}(\xi_k)}\sum_{k=1}^n\int_{|x+\mathbb{E}(\xi_k)|\ge \tau\sqrt{\sum_{k=1}^n \mathbb{D}(\xi_k)}}(x-\mathbb{E}(\xi_k))^2\mathrm{d}F_k(x)\to 0, \]

就有

\[\frac{\sum_{k=1}^n(\xi_k-\mathbb{E}(\xi_k))}{\sqrt{\sum_{k=1}^n \mathbb{D}(\xi_k)}}\stackrel{d}\to N(0,1). \]

這說明自然界中微小隨機項的累積效應普遍服從中心極限定理。

另外，正態分布的信息完全由兩個參數所決定：期望和方差，即前兩階矩。因此，如果我們假定總體是服從正態分布的，就只需要對其兩個參數作估計，這給問題的討論帶來方便。最后就是正態分布在實用上的意義了，兩個獨立正態分布的和、差甚至乘積都是正態分布，這在實用上也很方便，所以許多時候即使總體不服從正態分布，也近似認為服從正態分布。

Part 2：正態分布均值估計

既然正態分布完全由兩個參數所決定，那么只要知道出這兩個參數的值（或者范圍），就能確定總體的全部信息。然而，在實際生活中要獲得絕對正確的正態分布參數是不可能的，因為生活中的總體情況總是未知，要認識總體，我們只能從總體中抽取一系列樣本，再通過樣本性質來估計總體。

最簡單的情況是簡單隨機抽樣，這時候每一個樣本都和總體具有相同的分布函數或密度函數。具體對於正態分布來說，\(X\sim N(\mu,\sigma^2)\)，如果我們抽取了\(n\)個簡單隨機樣本\((X_1,X_2,\cdots,X_n)\)，則\(X_1,\cdots,X_n\)之間實際上相互獨立，且\(\forall i,X_i\sim N(\mu,\sigma^2)\)。盡管\(\mu\)和\(\sigma^2\)我們未知，但是我們知道一點——它們一定是不會變化的常數，這樣，我們能夠獲得獨立且與總體分布相同的樣本，通過觀測樣本構造統計量來估計總體。這種將統計量的觀測值作為參數估計的估計方式，稱為點估計。

對於總體均值，很自然的一點是用樣本均值作為總體均值的估計。似乎沒有理由不這么做，但這么做有什么依據嗎？我們知道，觀測樣本具有兩重性，所以統計量也具有兩重性。要研究用樣本均值作為總體均值估計的合理性，必須觀察樣本均值作為隨機變量時的分布。

正態分布具有可加性，這指的是對於相互獨立的正態分布，它們的和作為一個隨機變量仍然服從正態分布，且均值和方差都是各分量的直接加和。有了這一點，我們就可以研究樣本均值的分布了。

\[T_n\xlongequal{def}\sum_{j=1}^nX_j \\ \bar X=\frac{1}{n}\sum_{j=1}^n X_j=\frac{1}{n}T_n, \]

由於正態分布服從可加性，因此有

\[T_n\sim N(n\mu,n \sigma^2). \]

另外，由於正態分布的數乘依然是正態分布，且均值相當於乘上常數，方差相當於乘上常數的平方，所以

\[\bar X=\frac{1}{n}T_n\sim N\left(\mu,\frac{\sigma^2}{n}\right). \]

直觀上來看，樣本均值與總體具有相同的均值，但是方差變成了原來的\(n\)分之一。眾所周知，方差代表隨機變量取值的離散情況，由切比雪夫不等式有\(\forall\varepsilon>0\)，

\[\mathbb{P}(|\bar X-\mu|>\varepsilon)\le \frac{\mathbb{D}(\bar X)}{\varepsilon^2}=\frac{\sigma^2}{n\varepsilon^2}. \]

這個式子表明，只要你的\(n\)無限增大，\(\bar X\)和\(\mu\)之間的差距就可以無限縮小，這個性質稱為弱相合性。另外，樣本均值的期望和要估計的參數\(\mu\)一致，這個性質稱為無偏性。由於樣本均值估計總體均值時具有無偏性、弱相合性等優點，所以我們很難不把樣本均值當作總體均值的點估計。

下面是用R語言從\(N(5,25)\)中獲得100次樣本均值的程序，每一個樣本均值是由10000個相互獨立的樣本構成的。

rm(list = ls())  # 清空內存
barxlst <- c()
for (i in 1:100){
  barxlst[i] <- mean(rnorm(10000, 5, 5))
}
split.screen(c(1, 2))
screen(1)
plot(barxlst)  # 繪制散點圖
screen(2)
hist(barxlst)  # 繪制直方圖

從圖上可以看到，樣本均值的集中區間幾乎都在\(4.9\sim 5.1\)之間，讀者可以自行利用正態分布的\(3\sigma\)性質驗證這一點。

Part 3：正態分布方差估計

接下來就輪到第二個參數\(\sigma^2\)了，大家也很容易想到用樣本方差

\[S^2=\frac{1}{n-1}\sum_{j=1}^n(X_j-\bar X)^2 \]

作為總體方差的估計，但隨之而來的就有兩個問題：

1. 為什么要用\(S^2\)這個估計量？
2. 為什么\(S^2\)的分母是\(n-1\)而不是\(n\)？

這里我們依然要探究\(S^2\)的分布，但在此之前，先探究一下\(S^2\)的均值，這比探究\(S^2\)的分布要更容易。為此，可以作如下的變形：

\[\begin{aligned} S^2&=\frac{1}{n-1}\sum_{j=1}^n(X_j-\bar X)^2\\ &=\frac{1}{n-1}\sum_{j=1}^n[X_j-\mu-(\bar X-\mu)]^2\\ &=\frac{1}{n-1}\left[\sum_{j=1}^n(X_j-\mu)^2+n(\bar X-\mu)^2-2\sum_{j=1}^n(X_j-\mu)(\bar X-\mu) \right]. \end{aligned} \]

對於每一個\(X_j\)，有\((X_j-\mu)\sim N(0,\sigma^2)\)，所以易得\(\mathbb{E}(X_j-\mu)^2=\sigma^2\)；由於\((\bar X-\mu)\sim N(0,\sigma^2/n)\)，所以\(\mathbb{E}(\bar X-\mu)^2=\sigma^2/n\)。最后一部分，有

\[\begin{aligned} &\quad \sum_{j=1}^n(X_j-\mu)(\bar X-\mu)\\ &=(\bar X-\mu)\left[\sum_{j=1}^nX_j-n\mu\right]\\ &=n(\bar X-\mu)^2, \end{aligned} \]

所以

\[\begin{aligned} \mathbb{E}(S^2)&=\frac{1}{n-1}\left[\sum_{j=1}^n\mathbb{E}(X_j-\mu)^2-n\mathbb{E}(\bar X-\mu)^2 \right]\\ &=\frac{n}{n-1}\sigma^2-\frac{n}{n-1}\cdot\frac{\sigma^2}{n}\\ &=\sigma^2. \end{aligned} \]

這說明\(S^2\)在估計\(\sigma^2\)上是無偏的，這樣我們就解決了第二個問題：為什么\(S^2\)的分母是\(n-1\)而不是\(n\)，這是因為\(n-1\)作分母可以讓統計量具有無偏性。

下面就是為什么要用\(S^2\)估計\(\sigma^2\)的問題，照理說具有無偏性的估計量可以有這么多，為什么選擇了樣本方差呢？回想總體均值的估計，事實上如果我們只想獲得一個無偏估計，使用\(X_1\)就夠了（顯然有\(\mathbb{E}(X_1)=\mu\)），但用\(X_1\)估計不具有相合性，也就是說不管你樣本容量多大，這個統計量都不向着真值的方向靠近，這顯然不是我們想要的效果。使用\(S^2\)是否也具有一樣的相合性，是我們需要驗證的問題。

注意下面的證明步驟十分重要，請大家務必將它記下來。

對於樣本\(X_1,\cdots,X_n\)，使用施密特正交化構造一個如下的正交陣（不會不知道什么叫正交陣了吧）：

\[A=\begin{bmatrix} \frac{1}{\sqrt{n}} & \frac{1}{\sqrt{n}} & \cdots & \frac{1}{\sqrt{n}} \\ a_{21} & a_{22} & \cdots & a_{2n} \\ \vdots & \vdots & & \vdots \\ a_{n1} & a_{n2} & \cdots & a_{nn} \end{bmatrix}, \]

這個正交陣是一定可以構造出來的，如果你覺得很玄，下面是一個構造方式：

\[A=\begin{bmatrix} \frac{1}{\sqrt{n}} & \frac{1}{\sqrt{n}} & \frac{1}{\sqrt{n}} & \dots & \frac{1}{\sqrt{n}} \\ \frac{1}{\sqrt{2\cdot 1}} & \frac{-1}{\sqrt{2\cdot 1}} & 0 & \cdots & 0 \\ \frac{1}{\sqrt{3\cdot 2}} & \frac{1}{\sqrt{3\cdot 2}} & \frac{-2}{\sqrt{3\cdot 2}} & \cdots & 0 \\ \vdots & \vdots & \vdots & & \vdots \\ \frac{1}{\sqrt{n(n-1)}} & \frac{1}{\sqrt{n(n-1)}} & \frac{1}{\sqrt{n(n-1)}} & \cdots & \frac{-(n-1)}{\sqrt{n(n-1)}} \end{bmatrix}. \]

令\(\boldsymbol{X}=(X_1,\cdots,X_n)'\)，

\[\boldsymbol{Y}=(Y_1,\cdots,Y_n)'\xlongequal{def}AX, \]

則

\[Y_1=\sum_{j=1}^n\frac{1}{\sqrt{n}}X_j=\sqrt{n}\bar X. \]

由正交變換保持向量長度不變，得到

\[X_1^2+\cdots+X_n^2=Y_1^2+\cdots+Y_n^2, \]

所以

\[(n-1)S^2=\sum_{j=1}^n(X_j-\bar X)^2=\sum_{j=1}^nX_j^2-n\bar X^2=\sum_{j=1}^n Y_j^2-Y_1^2=\sum_{j=2}^nY_j^2. \]

接下來要證明一個很神奇的結論：\(Y_2,\cdots,Y_n\)獨立同分布於\(N(0,\sigma^2)\)。首先由正態分布的線性組合仍然是正態分布這一性質，知道\(Y_2,\cdots,Y_n\)都服從正態分布，而它們的均值，對\(i=2,\cdots,n\)，有

\[\mu_i=\sum_{j=1}^na_{ij}\mathbb{E}(X_j)=\sum_{j=1}^na_{ij}\mu=\sqrt{n}\mu\sum_{j=1}^n a_{ij}\frac{1}{\sqrt{n}}=\sqrt{n}\mu\sum_{j=1}^n a_{ij}a_{1j}=0. \]

這里用到了正交矩陣第\(j\)行和第\(1\)行的正交性。它們的方差，對\(i=2,\cdots,n\)，有

\[\sigma_i^2=\sum_{j=1}^na_{ij}^2\mathbb{D}(X_j)=\sigma^2\sum_{j=1}^na_{ij}^2=\sigma^2, \]

這里用到了正交矩陣每一行平方和為1的性質。接下來還要證明\(\forall j\ne k\)，\(Y_jY_k\)相互獨立（對於正態分布，即不相關），有

\[\mathrm{Cov}(Y_j,Y_k)=\sum_{i=1}^n a_{ij}a_{ik}\sigma^2=\sigma^2\sum_{i=1}^na_{ij}a_{ik}=0. \]

這就說明\(Y_1,Y_2,\cdots,Y_n\)相互獨立，且\(Y_2,\cdots,Y_n\)獨立同分布於\(N(0,\sigma^2)\)，所以它們的平方和為

\[(n-1)S^2=\sum_{j=2}^nY_j^2,\\ S^2=\frac{1}{n-1}\sum_{j=2}^n Y_j^2 \]

這里

\[\mathbb{E}(Y_j^2)=\sigma^2,\quad \mathbb{D}(Y_j^2)=\mathbb{E}(X_j^4)-\sigma^4=2\sigma^4,\\ \begin{aligned} \mathbb{E}(X_j^4)&=\int_{-\infty}^\infty x^4\frac{1}{\sqrt{2\pi\sigma^2}}e^{-\frac{x^2}{2\sigma^2}}\mathrm{d}x\\ &=\frac{1}{\sqrt{2\pi\sigma^2}}\int_{-\infty}^\infty x^4\frac{-x}{\sigma^2}\frac{\sigma^2}{-x}\exp\left(-\frac{x^2}{2\sigma^2} \right)\mathrm{d}x\\ &=\frac{1}{\sqrt{2\pi\sigma^2}}(-\sigma^2)\int_{-\infty}^\infty x^3\mathrm{d}(e^{-\frac{x^2}{2\sigma^2}})\\ &=-\frac{\sigma}{\sqrt{2\pi}}\left(x^3e^{-\frac{x^2}{2\sigma^2}}\bigg|_{-\infty}^\infty -\int_{-\infty}^{\infty}3x^2e^{-\frac{x^2}{2\sigma^2}}\mathrm{d}x \right)\\ &=\frac{3\sigma^2}{\sqrt{2\pi\sigma^2}}\int_{-\infty}^{\infty}x^2e^{-\frac{x^2}{2\sigma^2}}\mathrm{d}x\\ &=3\sigma^2\mathbb{E}(X_j^2)\\ &=3\sigma^4. \end{aligned} \]

所以

\[\mathbb{D}(S^2)=\frac{1}{n-1}\cdot\mathbb{D}(Y_1^2)=\frac{2\sigma^4}{n-1}, \\ \mathbb{P}(|S^2-\sigma^2|>\varepsilon)\le \frac{\mathbb{D}(S^2)}{\varepsilon^2}=\frac{2\sigma^4}{(n-1)\varepsilon^2}. \]

這就證明了\(S^2\)估計\(\sigma^2\)的弱相合性。同時，由於\(\bar X\)只是\(Y_1\)的函數，\(S^2\)只是\(Y_2,\cdots,Y_n\)的函數，由\(Y_1\)與\(Y_2,\cdots,Y_n\)的相互獨立性，我們還能得到一個重要結論：\(\bar X\)與\(S^2\)相互獨立。

以上的論證說明，\(\bar X\)作為\(\mu\)的估計是無偏且相合的，\(S^2\)作為\(\sigma^2\)的估計是無偏且相合的，並且\(\bar X\)與\(S^2\)相互獨立。

R語言中，var()函數能獲得樣本方差，對\(N(5,25)\)的10000個樣本、1000000個樣本進行100次模擬，注意觀察坐標軸。

rm(list = ls())
varlst1 <- c()
varlst2 <- c()
for (i in 1:100){
  varlst1[i] <- var(rnorm(10000, 5, 5))
  varlst2[i] <- var(rnorm(1000000, 5, 5))
}

split.screen(c(2, 2))

screen(1)
plot(varlst1)
title("10000個樣本：散點圖")

screen(2)
hist(varlst1, main = "10000個樣本：直方圖")

screen(3)
plot(varlst2)
title("1000000個樣本：散點圖")

screen(4)
hist(varlst2, main = "1000000個樣本：直方圖")
# dev.off()

Part 4：卡方分布

以上對探索\(S^2\)的分布終究只是從它的一階矩、二階矩上證明了它用於刻畫\(\sigma^2\)的優良性質，有沒有辦法能夠得到\(S^2\)分布的詳細信息呢？我們來觀察\(S^2\)的表達式，注意到

\[(n-1)S^2=\sum_{j=2}^n Y_j^2=\sigma^2\sum_{j=2}^n\left(\frac{Y_j}{\sigma} \right)^2,\\ \frac{(n-1)S^2}{\sigma^2}=\sum_{j=2}^n\left(\frac{Y_j}{\sigma} \right)^2. \]

這樣變換的意義在於，右邊變成了\(n-1\)個獨立同分布的標准正態隨機變量的平方和，在數理統計中，我們會經常遇到類似於這樣的分布，因此將其定義為卡方分布。以后我們會了解到，正態分布有三大常用的衍生分布，今天只介紹第一種。

定義：設\(X_1,\cdots,X_n\)獨立同分布於\(N(0,1)\)，則稱

\[\xi=\sum_{j=1}^n X_j^2\sim \chi^2(n). \]

這里\(n\)稱為\(\chi^2\)分布的自由度。

從我們剛才對\(S^2\)的討論，容易知道如果\(\xi\sim \chi^2(n)\)，則

\[\mathbb{E}(\xi)=n,\quad \mathbb{D}(\xi)=2n. \]

進一步，我們還可以求出\(\xi\)的密度函數為

\[f_\xi(x)=\frac{1}{2^{\frac{n}2}\Gamma(\frac{n}{2})}x^{\frac{n}{2}-1}e^{-\frac{x}{2}}I_{(x>0)}. \]

這里\(I_A\)是示性函數，即\(A\)發生時取1不發生時取0。在數理統計中，選擇示性函數而不是分段函數的形式表示密度函數會給問題的討論帶來極大的方便。目前，記憶卡方分布的密度還挺有難度的，但我們可以暫且跳過它。

由於\(X_1,\cdots,X_n\)獨立同分布服從於\(N(0,1)\)，所以其聯合密度是

\[f(x_1,\cdots,x_n)=\left(\frac{1}{\sqrt{2\pi}} \right)^n\exp\left\{-\frac{1}{2}\sum_{j=1}^nx_j^2 \right\}, \]

設\(\xi=\sum_{j=1}^n X_j^2\)的分布函數為\(F_n(x)\)，則

\[F_n(x)=\left(\frac{1}{\sqrt{2\pi}} \right)^2\int\cdots\int_{\sum_{j=1}^nx_j^2<x}\exp\left\{-\frac{1}{2}\sum_{j=1}^nx_j^2 \right\}\mathrm{d}x_1\cdots\mathrm{d}x_n, \]

作以下球坐標變換：

\[\left\{\begin{array}l x_1=\rho\cos(\theta_1)\cos(\theta_2)\cdots\cos(\theta_{n-2})\cos(\theta_{n-1}), \\ x_2=\rho\cos(\theta_1)\cos(\theta_2)\cdots\cos(\theta_{n-2})\sin(\theta_{n-1}), \\ \vdots \\ x_n=\rho\sin(\theta_1), \end{array}\right. \]

該變換的Jacobi行列式絕對值為

\[|J|=\left|\frac{\partial(x_1,x_2,\cdots,x_n)}{\partial(\rho,\theta_1,\cdots,\theta_{n-1})} \right|=\rho^{n-1}g(\boldsymbol{\theta}). \]

這里\(\rho\le \sqrt{x}\)，\(g(\boldsymbol{\theta})\)是某個關於\(\boldsymbol{\theta}=(\theta_1,\cdots,\theta_{n-1})\)的函數，所以

\[\begin{aligned} F_n(x)&=\left(\frac{1}{\sqrt{2\pi}} \right)^2\int_{0}^{\sqrt{x}}\rho^{n-1}e^{-\frac{\rho^2}{2}}\mathrm{d}\rho\cdot\int_{\Theta} g(\boldsymbol{\theta})\mathrm{d}\boldsymbol{\theta}\\ &=C\int_{0}^{\sqrt{x}}\rho^{n-1}e^{-\frac{\rho^2} {2}}\mathrm{d}\rho\\ &\xlongequal{t:=\rho^2}C\int_{0}^{x}t^{\frac{n}{2}-1}e^{-\frac{t}{2}}\mathrm{d}t. \end{aligned} \]

這里\(C\)是某個常數，后面的部分為關於\(t\)的核，故顯然\(\xi\)的密度為

\[f_n(x)=Cx^{\frac{n}{2}-1}e^{-\frac{x}{2}}I_{x>0}, \]

后面的部分稱為密度函數的核，只要確定常數\(C\)即可，經過積分得到

\[C=\frac{1}{2^{\frac{n}{2}}\Gamma(\frac{n}{2})}. \]

最后指出\(\chi^2\)具有可加性，這指的是對於獨立的\(X\sim \chi^2(n_1)\)，\(Y\sim \chi^2(n_2)\)，有

\[X+Y\sim \chi^2(n_1+n_2), \]

只要把\(X\)看成\(n_1\)個獨立標准正態變量的平方和，\(Y\)看成\(n_2\)個獨立標准正態變量的平方和，它們的和就是\(n_1+n_2\)個獨立標准正態變量的平方和，故服從\(\chi^2(n_1+n_2)\)。

結合卡方分布，我們可以對今天得出的結論作一個小小的總結了。設\(X\sim N(\mu,\sigma^2)\)中抽取的簡單隨機樣本導出的樣本均值為\(\bar X\)，樣本方差為\(S^2\)，則

1. \(\bar X\)的分布：

   \[\bar X\sim N\left(\mu,\frac{\sigma^2}{n} \right). \]

2. \(S^2\)的分布：

   \[\frac{(n-1)S^2}{\sigma^2}\sim \chi^2(n-1). \]

3. \(\bar X\)與\(S^2\)相互獨立。

---

最后，有一個問題需要大家思考：正態參數的總體均值\(\mu\)和總體方差\(\sigma^2\)的無偏、相合估計量只有\(\bar X\)和\(S^2\)嗎？如果不是，為什么我們會選擇\(\bar X\)和\(S^2\)來估計總體參數呢？舉個最簡單的例子，如果我們取出了\(n\)個樣本，那么用前\(n/2\)個（如果\(n\)是偶數）樣本來計算樣本均值和樣本方差，一樣是相合且無偏的。問題在於，為什么要用全部\(n\)個樣本來計算，而不是只用部分的樣本。

之后，我們會對更多的分布族討論參數估計問題。