上篇文章中,我們探討了區間估計的相關基本概念,也提出了Neyman置信區間,今天我們將聚焦於如何尋找置信區間的問題上,並對最常用的總體:正態總體給出一些置信區間的找法。為了方便起見,以下我們都讓置信水平為\(1-\alpha\)。
由於本系列為我獨自完成的,缺少審閱,如果有任何錯誤,歡迎在評論區中指出,謝謝!
Part 1:樞軸量法
樞軸變量法是基於點估計量的。我們知道,統計量是樣本的函數,這意味着統計量中不能含有未知參數,而參數的點估計量是用統計量的觀測值作為待估參數的估計值,其分布一定含有待估參數,樞軸量法的思想就是,通過一定的變換,讓點估計的函數的分布不含待估參數,進而基於分布來構造區間估計。
舉一個簡單的例子,對於正態總體\(N(\mu,4)\),顯然\(\bar X\sim N(\mu,4/n)\),這里\(\bar X\)的分布含有未知參數\(\mu\)。構造其樞軸量,就是找到一個函數變換,使得新的隨機變量分布不含未知參數。注意,這里用了隨機變量這個詞而不是統計量,意味着樞軸量不是統計量,即不能由樣本觀測值計算出,這是因為雖然樞軸量的分布不含未知參數,但是樞軸量的表現形式含有未知參數。顯然,這里
這樣,\(\bar X-\mu\)的分布已知,自然容易找到一個常數區間\([c,d]\),使得這個區間有\(1-\alpha\)的概率包含\(\bar X-\mu\)的觀測值,雖然此時我們不知道區間的端點是多少,但至少知道端點可以是固定的數\(c,d\)。對樞軸量使用不等式變換,即\(\bar X-\mu\in[c,d]\Rightarrow \mu\in[\bar X-d,\bar X-c]\),得到置信水平為\(1-\alpha\)的置信區間。這就是樞軸量法的操作步驟。
不同分布族的參數對於總體的意義是不同的。像正態分布\(N(\mu,\sigma^2)\)的均值\(\mu\),均勻分布\(U(a,a+r)\)的起點\(a\)這種參數主要影響觀測值的大小,可以直接通過\(X-\mu\),\(X-a\)的手段消除,這種參數稱為位置參數;正態分布\(N(\mu,\sigma^2)\)的標准差\(\sigma\),指數分布\(E(\lambda)\)的速率\(\lambda\)這種參數主要影響觀測值的離散程度,可以通過\(X/\sigma\),\(\lambda X\)之類的手段消除,這種參數稱為尺度參數。面對不同的參數,往往有不同的方式構造樞軸量,應當結合其特點選擇構造方式。
接下來,我們將樞軸量法運用到一些具體實例中,體會樞軸量法的實際應用。
Part 2:分位數
在開始樞軸量法的實際應用前,先介紹分位數這一概念,這為我們確定置信區間提供了有力武器。現在,我們先給出分位數的定義(如果不特別說明,都指總體分位數而非樣本分位數)。分位數可以分為下分位數和上分位數,一般稱下分位數為分位數,但我們更常用的是上分位數。
對於某個分布\(F\),\(X\sim F\),\(F\)的(下)\(\alpha\)分位數指的是這樣一個點\(x_\alpha\),滿足
如果\(F\)有反函數\(F^{-1}(x)\),則有\(x_\alpha=F^{-1}(\alpha)\)。
對於某個分布\(F\),\(X\sim F\),\(F\)的上\(\alpha\)分位數指的是這樣一個點\(y_\alpha\),滿足
換言之,如果分布函數\(F\)存在反函數\(F^{-1}(x)\),則\(F\)的上\(\alpha\)分位數是\(y_\alpha=F^{-1}(1-\alpha)\)。

也許結合圖形會更好理解。對於上面的圖形,如果這是一個密度函數,黑色部分的面積為\(0.05\),那么紅色的點就是該分布的上\(0.05\)分位數,同時也是其下\(0.95\)分位數。
顯然,對於一個完全確定的分布,其各分位數都是完全確定的,是常數。分布多種多樣,不可能對所有分布都詳細地討論其分位數,但是對一些常用分布,其各分位數的數值已經被制成表,這包括標准正態分布,與正態分布的三大衍生分布——\(\chi^2\)分布、\(t\)分布、\(F\)分布。
今后,我們將使用\(u_\alpha\)、\(\chi^2_\alpha(n)\)、\(t_\alpha(n)\)、\(F_{\alpha}(m,n)\),分別代表標准正態分布、自由度為\(n\)的\(\chi^2\)分布、自由度為\(n\)的\(t\)分布、自由度為\(m,n\)的\(F\)分布的上\(\alpha\)分位數,給定分布類型、分位、自由度后,這些值都可以通過查閱分位數表得到。
Part 3:單正態總體參數區間估計
正態分布參數的區間估計,主要是對\(\mu\)和\(\sigma^2\)給出區間估計。其中,對\(\mu\)的估計又要分成\(\sigma^2\)已知和未知兩種情況。不過,大家已經知道,與參數\(\mu\)關聯最緊密的點估計是\(\bar X\),與參數\(\sigma^2\)關聯最緊密的點估計是\(S^2\),因此樞軸量法也一定圍繞着它們作文章。
事實上,只要是均值,基本都可以轉化為一個正態分布;只要是方差,基本上都可以轉化為一個卡方分布。然后,根據\(t\)分布和\(F\)分布的構造方式,就能構造出服從這四種分布的樞軸量來。
1、\(\sigma^2\)已知時\(\mu\)的區間估計
當\(\sigma^2\)已知時,
樞軸量顯然是\(\bar X-\mu\sim N(0,\sigma^2/n)\),已經由此確定了一個已知分布,但是為了寫出區間估計的具體表達,我們還應當將這個已知的正態分布,轉化為標准正態分布。為此,我們選擇樞軸量為
找一個以\(1-\alpha\)概率涵蓋標准正態分布觀測的區間,自然會找到
即
進行恆等變換,就得到了\(\mu\)的\(1-\alpha\)置信區間為
2、\(\sigma^2\)未知時\(\mu\)的區間估計
當\(\sigma^2\)未知時,我們不能順利地將\(\bar X\)標准化,但應該容易想到,用\(S^2\)替代\(\sigma^2\)能起到類似的效果,也就是構造這樣的變量:\(\sqrt{n}(\bar X-\mu)/S\)。
神奇的是,\(\frac{\sqrt{n}(\bar X-\mu)}{S}\)可以被改寫成\(\chi^2\)分布的形式:
因此,自然可以找到一個區間\(\frac{\sqrt{n}(\bar X-\mu)}{S}\in[-t_{\alpha/2}(n-1),t_{\alpha/2}(n-1)]\),因此得到\(\mu\)的\(1-\alpha\)置信區間為
用R語言進行模擬,現在我們每次生成10個服從\(N(5,9)\)的簡單隨機樣本,這一組樣本可以生成一個對應的\(95\%\)置信區間(假設\(\sigma^2=9\)是未知的)。進行10000次試驗,生成10000個這樣的區間,觀察這個區間會涵蓋待估參數\(\mu=5\)多少次。
rm(list = ls()) # 清空工作區
dev.off() # 清空圖形區
n <- 10
t <- qt(0.975, df=n-1) # t的上0.025分位數
upper.bound <- c()
lower.bound <- c()
for (j in 1:10000){
xlst <- rnorm(10, 5, 3)
M <- mean(xlst)
S <- sqrt(var(xlst))
upper.bound[j] <- M+S*t/sqrt(n-1)
lower.bound[j] <- M-S*t/sqrt(n-1)
}
plot(1:10000, rep(5, 10000), main = "Mean of Norm", ylim = c(-10, 20), cex=0.1)
points(1:10000, upper.bound, col='red', cex=0.1)
points(1:10000, lower.bound, col='blue', cex=0.1)
sum(upper.bound < 5 | lower.bound > 5) # 輸出區間不涵蓋待估參數的次數
輸出結果為,10000次中,有384次沒有涵蓋待估參數,圖示如下:

3、\(\mu\)未知時\(\sigma^2\)的區間估計
我們在給出\(S^2\)的分布信息時,實際上已經構建出了樞軸量:
所以我們可以給出相應的\(1-\alpha\)概率區間,不過要注意的是\(\chi^2\)分布不像正態分布、\(t\)分布一樣,是對稱分布,所以概率區間為\(\frac{(n-1)S^2}{\sigma^2}\in[\chi^2_{1-\alpha/2}(n-1),\chi^2_{\alpha/2}(n-1)]\),得到\(\sigma^2\)的置信區間為
使得概率為\(1-\alpha\)的區間理論上有無窮多個,它們有相同的置信度,本該選擇其中精度最高(即長度最短)的,但是這樣會很麻煩,為方便起見,選取\(\alpha/2\)分位數和\(1-\alpha/2\)分位數構造置信區間,這樣的置信區間稱為等尾置信區間。
4、\(\mu\)已知時\(\sigma^2\)的區間估計
看起來,\(S^2\)這個樞軸量不考慮\(\mu\)是否已知,但是樞軸量法是依賴於點估計的,自然我們首先會考慮的是充分統計量,而\(S^2\)在\(\mu\)已知時,不是\(\sigma^2\)的充分統計量。為此,我們要基於\(\sigma^2\)的充分統計量構建樞軸量。令
它是\(\sigma^2\)的充分統計量,且我們可以將其轉化為\(\chi^2\)分布:
這是由於各個\(X_j\)相互獨立服從\(N(\mu,\sigma^2)\),於是可以構造出概率區間\(\frac{Q}{\sigma^2}\in[\chi^2_{1-\alpha/2}(n),\chi^2_{\alpha/2}(n)]\),因此\(\sigma^2\)的置信區間是
Part 4:雙正態總體參數區間估計
雙總體\(X\sim N(\mu_1,\sigma_1^2)\)、\(Y\sim N(\mu_2,\sigma_2^2)\)的參數估計,主要指的是均值差和方差比,但是要分情況討論,這是因為如果不加任何限制,它們的置信區間理論還不是很完善。
以下記\(X_1,\cdots,X_m\)是從\(X\)中抽取的樣本,相應的樣本均值和樣本方差是\(\bar X\)和\(S_x^2\);\(Y_1,\cdots,Y_n\)是從\(Y\)中抽取的樣本,相應的樣本均值和樣本方差是\(\bar Y\)和\(S_y^2\)。
1、均值差\(\mu_2-\mu_1\)的區間估計
第一種情況是\(\sigma_1^2,\sigma_2^2\)未知但\(m=n\)時,即樣本容量相等時。一種簡單的做法是構造\(Z_j=Y_j-X_j\),即讓成對數據相減,這樣\(Z_1,\cdots,Z_n\)便獨立同分布於\(N(\mu_2-\mu_1,\sigma_1^2+\sigma_2^2)\),相當於從單總體中抽樣,並對總體均值作估計。顯然,樞軸量應該是
所以置信區間是
不過要注意的是,構造成對數據時,實際上丟失了部分信息,這是因為\(Z_1,\cdots,Z_n\)並非\(\mu_2-\mu_1\)的充分統計量。
總而言之,當\(m=n\)時,我們將其轉化為成對數據,本質上還是單總體參數估計。
第二種情況是\(\sigma_1^2,\sigma_2^2\)已知時,此時用\(\bar Y-\bar X\)來最大程度地估計\(\mu_2-\mu_1\),則有
於是樞軸量是
\(\mu_2-\mu_1\)的置信區間為
第三種情況稍微復雜一些,也是最常處理的一種情況,即\(\sigma_1^2=\sigma_2^2=\sigma^2\)但未知時。此時仍有
為了解決\(\sigma^2\)未知的問題,也要用適當的樣本方差來替代,但此時如此構造:
可以發現
於是用\(Q/(n+m-2)\)代替\(\sigma^2\)最合適不過了,因此
即\(\mu_2-\mu_1\)的置信區間為
現進行模擬,每次實驗中,從\(N(3,9)\)中抽取5個樣本,從\(N(5,25)\)中抽取8個樣本,對均值差進行估計。進行10000次重復試驗。
rm(list = ls())
dev.off()
m <- 5
n <- 8
t <- qt(0.975, m+n-2)
upper.bound <- c()
lower.bound <- c()
for (j in 1:10000){
xlst <- rnorm(m, 3, 3)
ylst <- rnorm(n, 5, 5)
My <- mean(ylst)
Mx <- mean(xlst)
Q <- (m-1)*var(xlst) + (n-1)*var(ylst)
upper.bound[j] <- My-Mx+t*sqrt( (Q*(m+n)) / ((m+n-2)*m*n) )
lower.bound[j] <- My-Mx-t*sqrt( (Q*(m+n)) / ((m+n-2)*m*n) )
}
plot(1:10000, rep(2,10000), xlab = "Experiment", ylab = "CI", main = "CI of the mean difference", cex=0.1, ylim = c(-13, 17))
points(1:10000, upper.bound, col = "red", cex = 0.1)
points(1:10000, lower.bound, col = "blue", cex = 0.1)
sum(upper.bound < 2 | lower.bound > 2) # 輸出結果為302

當我們將樣本容量提升至\(m=200\),\(n=100\)時,出現了929個不包含真值的區間。看起來區間估計是否包含真值,與\(m,n\)的相對大小有很大的關系。

當\(\sigma_1^2\ne \sigma_2^2\)且未知時,均值差的參數估計稱為Behrens-Fisher問題,比較復雜,這里不加討論。
2、方差比\(\sigma_2^2/\sigma_1^2\)的區間估計
事實上方差比比較好討論,因為方差的估計量一定是自己總體中所抽取的樣本方差,故服從某個\(\chi^2\)分布,作比后應當能得到\(F\)分布,這里簡要地討論兩種情況。
首先是\(\mu_1,\mu_2\)都未知的情況,此時自然\(S_x^2,S_y^2\)會被用作估計量,由於
所以
故應有\(\frac{S_x^2\cdot\sigma_2^2}{S_y^2\cdot\sigma_1^2}\in[F_{1-\alpha/2}(m-1,n-1),F_{\alpha/2}(m-1,n-1)]\),所以置信區間為
進行模擬,每次實驗中,從\(N(3,9)\)中抽取9個樣本,從\(N(5,25)\)中抽取25個樣本,對方差比進行估計。進行10000次重復試驗。
rm(list = ls())
dev.off()
m <- 9
n <- 25
f1 <- qf(0.025, m-1, n-1)
f2 <- qf(0.975, m-1, n-1)
upper.bound <- c()
lower.bound <- c()
for (j in 1:10000){
xlst <- rnorm(m, 3, 3)
ylst <- rnorm(n, 5, 5)
Sx2 <- var(xlst)
Sy2 <- var(ylst)
lower.bound[j] <- Sy2/Sx2 * f1
upper.bound[j] <- Sy2/Sx2 * f2
}
plot(1:10000, rep(25/16,10000), xlab = "Experiment", ylab = "CI", main = "CI of the variance ratio", cex=0.1, ylim = c(0, 40))
points(1:10000, upper.bound, col = "red", cex = 0.1)
points(1:10000, lower.bound, col = "blue", cex = 0.1)
sum(upper.bound < 25/16 | lower.bound > 25/16)

然后是\(\mu_1,\mu_2\)均已知的情況,這時方差的估計量會是
顯然有
於是樞軸量為
即\(\sigma_2^2/\sigma_1^2\)的區間估計為
事實上,方差的估計,就是在\(\mu\)已知時使用對均值的離差平方和,在\(\mu\)未知時使用對樣本均值的離差平方和,注意分母與\(\chi^2\)分布的自由度即可。由此,讀者應該也能推斷出\(\mu_1\)已知\(\mu_2\)未知,或者\(\mu_1\)未知\(\mu_2\)已知時方差比的區間估計了。
本文給出了區間估計的求法——樞軸量法的實際應用案例,並對正態總體的參數作區間估計。下篇文章將着眼於非正態總體的區間估計,由於非正態總體不像正態總體一樣,擁有性質良好、分布透明的點估計可以使用,所以我們可能需要尋找一種更通用的方法來求區間估計。