有限群體之隨機抽樣與中央極限定理－劉應興的部落格

考慮有限群體 Y = {Y_1,...,Y_N}，此群體之一個大小為 n 的（簡單）隨機樣本 {X1,...,Xn} 是指從 Y 中取 n 個相異元素的 C(N,n) 種組合之一。此處「相異」指的是個體，意思是：從數值上可能 Y_i = Y_j for some i≠j，仍把相同數值的 Y_i, Y_j 看成是相異元素。或許更完整的數學描述是 P = {(i,Y_i): i=1,...,N}, 而樣本是

{Y_i t_i: t_i = 0 or 1, i = 1, ..., N, Σt_i = n}

其中隨機向量 (t_1, ..., t_N) 的機率分布是：

p(t_1,...,t_N) = 1/C(N,n) for t_i = 0 or 1, i = 1, ..., N, Σt_i = n

設 a_i 是 0 或 1，i = 1, ..., N，則特定 k 個 t_i = a_i 的機率，例如

P[t_i = a_i, i=1,...,k] = C(N-k, n-Σ_{i=1~k} a_i)/C(N,n)

例如：P[t_i = 1] = C(N-1,n-1)/C(N.n) = n/N, P[t_i=0] = 1-n/N，又設 i≠j，則

P[t_1 = 1, t_j = 1] = C(N-2,n-2)/C(N,n) = n(n-1)/[N(N-1)]

以上的意思是說：有限群體 P 的簡單隨機樣本是自 P 中整批抓取 n 個個體 {Yi(1), Yi(2), ..., Yi(n)} 形成的。但在實施上，簡單隨機樣本可能是自 P 一次抽取一個，抽出後不放回的方式抽取。也就是先自 P 中抽取 Yi(1)，而後自去除 Yi(1) 的群體中抽取 Yi(2)，再自去除 Yi(1), Yi(2) 的群體中抽取 Yi(3)，以此類推，至抽取 Yi(n) 結束。前者（整批抽取）相當於把 N 個個體隨機排序後抓取前 n 個，有 N! 種排序，但前 n 個及後 N-n 個個體各自順序 n!(N-n)! 種，不影響樣本內容，所以共有 C(N,n) 種不同的樣本組成。後者（不放回法逐次抽取）相當於 N 個個體隨機排序，留下第一個；剩下 N-1 個個體再隨機排序，再留下一個；剩下的又隨機排序，又留一個；以此類推，直至留下 n 個，剩下 N-n 個，則停止抽樣，留下的那 n 個個體即是我們的隨機樣本，故有 P(N,n) = N!/(N-n)! 種結果，但只論組成不論抽出順序則也是 C(N,n) = P(N,n)/n! 種不同的可能樣本，實際得到的樣本是 C(N,n) 種可能樣本之一。

用 t_i 標示 Yi (或 Y_i) 是否在樣本中，樣本平均數 Xbar = ( ΣXj)/n = (ΣY_i t_i)/n 之期望值為

E[Xbar] = (Σ Y_i E[t_i])/n = (Σ Y_i n/N)/n = ΣY_i/N = μ

而其變異數為

Var[Xbar] = (1/n^2){Σ_i Y_i^2 Var(t_i) + 2ΣΣ_{i<j} Y_i Y_j Cov(t_i,t_j)}

由 t_i 之邊際分布及 t_i, t_j 之聯合分布，

Var[t_i] = (n/N)(1-n/N), Cov(t_i,t_j) = -(n/N)(1-n/N)/(N-1)

定義 σ^2 = Σ(Y_i-μ)^2/N, S^2 = Nσ^2/(N-1)，則可得

Var[Xbar] = (S^2/n)(1-n/N) = (σ^2/n)[(N-n)/(N-1)]

與無限群體相比，S^2 相當於 σ^2，而 Var[Xbar] 多了 1-n/N 一個乘數，稱為：有限群校正數 (finite population correction factor)。有一個許多初學者常疑惑：究竟前述 fpc 是 1-n/N = (N-n)/N 或 (N-n)/(N-1)？由上式可知：問題在群體變異數是 S^2 或 σ^2？若用 S^2，則 fpc 是 (N-n)/N；若用 σ^2 則 fpc 是 (N-n)/(N-1)。再者，樣本變異數

s^2 = Σ(Xi - Xbar)^2/(n-1)

的期望值可證得是 S^2，因此若用樣本變異數估計群體變異數，則在估計 Var[Xbar] 時，fpc 用 (N-n)/N 較適當。當然，實務上 N 不會太小，fpc 的除數用 N 或 N-1 結果的差異幾乎可以忽略，前面的選擇只是就數學上的結果和統計估計的合理性而論。

樣本變異數 s^2 = Σ Xi^2/(n-1) - n(Xbar)^2/(n-1)，因此，

E[s^2] = E[Σ Y_i^2 t_i]/(n-1) - n(μ^2 + Var[Xbar])/(n-1)

將 E[t_i] = n/N 及前面 Var[Xbar] 結果代入，即可得 E[s^2] = S^2。

將 n/N 用 f 表示；又群體中心動差可定義成

μ_k = Σ((Yi - μ)^k/N

故 μ_2 = σ^2。則

E[(Xbar - μ)^2] = (1/n)(1-f)[N/(N-1)]μ_2

E[(Xbar - μ)^3] = (1/n^2)f(1-f)(1-2f){1-(N-4)/[(N-1)(N-2)]}μ_3

E[(Xbar - μ)^4]
= (3/n^2)(1-f)^2．
{([N(N^2-3N+3)]/[(N-1)(N-2)(N-3)])μ_2^2 - (N/[(N-2)(N-3)])μ_4}
+ (1/n^3)(1-f)(1-3f+f^2)．
{-12μ^4 - 12[(N-6)/(N-1)]μ^2μ_2 + 24((N-4)/[(N-1)(N-2)])μμ_3
- 3([N(N^2-5N+12)]/[(N-1)(N-2)(N-3)])μ_2^2
+ ((N^3+N^2-24N+72)/[(N-1)(N-2)(N-3)])μ_4}

故樣本平均數的抽樣分布具有偏態係數

(1/√n) (μ_3/μ_2^{3/2}) [f(1-2f)/√(1-f)] {(N^2-4N+6)/[(N-1)(N-2)]}

此結果可分四部分來看，首先是：當 n 增大時，Xbar 的偏態成 n^{-1/2} 等級縮小；其次，Xbar 的偏態和群體偏態係數 μ_3/μ_2^{3/2} 有直接關係；但抽樣比例 f = n/N 也影響 Xbar 的偏態高低，甚至抽樣比例太高，f > 1/2 時扭轉了偏斜方向；最後，群體大小當然也有些影響，不過除非 N 很小，其影響不大，最後一個因子亦可寫成 1 - (N-4)/[(N-1)(N-2)] 可能更容易看出其效果。

前面的樣本平均數第四階動差公式很複雜，我們把它分為兩部分，可以看出第一部分接近

3E[(Xbar - μ)^2] = (3/n^2)(1-f)^2[N/(N-1)]^2μ_2^2

實際上差距是

(3/n^2)(1-f)^2 ({[N(N^2-3)]/[(N-1)^2(N-2)(N-3)]}μ_2^2 - }N/[(N-2)(N-3)]}μ_4)

大約是 (3/n^2)(1-f)^2(μ_2^2/N-μ_4/N^2)，整體算 O(1/(n^2N) 等級。另外，第二部分是 O(1/n^3) 等級，相對比前一部分差距大。而 Xbar 的變異數是 O(1/n) 等級，所以，其峰度係數與常態峰 3 之差是 O(1/n) 等級，也就是說樣本數較大時，Xbar 的分布，偏態係數以 O(1/√n) 等級收斂到 0；峰度係數以 O(1/n) 等級收斂到常態峰。

假設有限群體 Y 實際上是自一無限群體 S 隨機抽取的一個樣本，其元素 Y_1, ..., Y_N 是 i.i.d. 的隨機變數的值，自 Y 中再以不放還法抽出 X_1, ..., X_n。如前述，樣本 X_1, ..., X_n 相當於把 Y_1, ..., Y_N 隨機排列之後取前 n 個，因此 X_1, ..., X_n 等價於從 S 抽出之隨機樣本。當群體 S 具二階動差時，依 i.i.d. 情形的中央極限定理，T_n = √n(Xbar - E[Y_1]) 向 N(0, Var(Y_1)) 做分布收斂，但這是在假設 Y_1, ..., Y_N 是自 S 抽出之隨機樣本，並且允許 Y_1, ..., Y_N 自由（隨機）變動時 T_n 的邊際分布。我們有個疑問：在給定 Y_1, ..., Y_N 的值之下，T_n 的條件分布是否也向某個 N(0,τ^2) 做分布收斂？由上述 Xbar 的偏態和峰度來看，這似乎是對的。但一個機率分布並不能由前四階動差唯一決定，因此前述偏態與峰度的觀察只能佐參，不能當做證明。

前面說過：自有限群體 Y 抽取的隨機樣本 X1,...,Xn，可以看成先把 Y_1,...,Y_N 做隨機排列，然後取前 n 個 Y 元素構成樣本。這在實作方面可以自 U(0,1) 無限群體抽取 U_1,...,U_N 這樣一組隨機樣本。令 U_i 對應 Y_i，R_i 為 U_i 的 rank（秩），即：將 U_i, i=1,...,N, 由小而大排序得順序統計量

U(1) < U(2) < . . . < U(N)

而後，若 U(j) = U_i, 則 R_i = j 是 U_i 的秩，故

t_i = 1 若且唯若 R_i ≦ n

即 Y_i 在樣本 {X1,...,Xn} 中若且唯若 R_i ≦ n。考慮另一個樣本 W = {W_1, W_2, ...，　W_m}，取樣規則如下：

Y_i 在 W 中若且唯若 U_i < f = n/N

這樣本的觀測值個數 m 是隨機的，服從樣本數 N 成功率 f 的二項分布。兩樣本 W 和 X = {X1,...Xn} 有什麼關係呢？

X 樣本是取 R_i ≦ n 的 Y_i 為樣本，樣本數 n 固定；

W 樣本是取 R_i ≦ m 的 Y_i 為樣本，樣本數 m 是 bin(N,f) 隨機變數。

所以，Ｗ樣本和 X 樣本相差在 m ≠ n。令

S_j = i 若且唯若 R_i = j

又令 X'j = Y_(S_j)。如 N = 5, 並且

U_3 < U_1 < U_5 < U_4 < U_2

則 S_j 為 (3, 1, 5, 4, 2), R_i 為 (2, 5, 1, 4, 3), 而

X' = (Y_3, Y_1, Y_5, Y_4, Y_2)

若 n = 2, m = 3, 則

X = {Y_3, Y_1}; W = {Y_3, Y_1, Y_5}

一般而言就是：

若 m < n, 則 W 樣本比 X 樣本少了 Y_{S_{m+1}}, ..., Y_{S_n};
若 m > n, 則 W 樣本比 X 樣本多了 Y_{S_{n+1}}, ..., Y_{S_m}.

令 T = Σ Xi, T' = Σ Ｗi, 則

若 m < n, 則 T - T' = Σ_{i=m+1 ～ n} X'_i = Σ_{i=m+1 ～ n} Y_{S_i};
若 m > n, 則 T' - T = Σ_{i=n+1 ～ m} X'_i = Σ_{i=n+1 ～ m} Y_{S_i}.

所以，當 m 已知且 m ≠ n 時，T 與 T' 的差相當於自 Y 群體抽出 |m-n| 個成員的隨機樣本總值，

E[(T-T')^2] = E[ E[(T-T')^2|m] ]
= E[Var(T-T'|m) + (E[T-T'|m])^2]
= E[|m-n|(1-|m-n|/N)[N/(N-1)]σ^2 + (m-n)^2μ^2]
≦ {√[Nf(1-f)] [N/(N-1)] - Nf(1-f)/(N-1)}σ^2 + Nf(1-f)μ^2

而原樣本值總和 T，即 n(Xbar) 的變異數是

Var(T) = n(1-f)[N/(N-1)]σ^2

應用公式 Var(T') = E[Var(T'|m)] + Var(E[T'|m]) 或直接由 T' = Σ Yi I_[Ui≦f] 各項相互獨立的事實，可得 T' 的變異數

Var(T') = n(1-f)(σ^2+μ^2)

因為考慮極限分布時，會把統計量做適當標準化，例如

T* = (T - nμ)/(√n σ); T'* = (T' - mμ)/(√n σ)

因此不妨假設 μ = ΣYi/N = 0。在 μ≠0 的情形，則是在計算 T 及 T' 時以 Yi-μ 取代 Yi。故

E[(T - T')^2]/Var(T) = 1/√[n(1-f)] - 1/N → 0

當 N, n 趨於無窮，而 n/N = f 保持固定或 f_N = n/N → f > 0 時。又，在不失一般性條件下假設 μ = 0, 則

Var(T)/Var(T') = N/(N-1) → 1

故 E[(T - T')^2]/Var(T') 同樣收斂至 0。事實上 Var(T) = [1+1/(N-1)]Var(T')，故

Var(T) - Var(T') = Var(T')/(N-1) → f(1-f)σ^2

可將原樣本總值的標準化量表示為

T/√Var(T) ＝ (T'/√Var(T'))(√Var(T')/√Var(T)) + (T-T')/√Var(T)

但 √Var(T')/√Var(T) → 1 且 (T-T')/√Var(T) 向 0 做均方收斂（因此也向 0 做機率收斂），故由 Slutsky's 定理，若 T'/√Var(T') 有極限分布 F，則 T/√Var(T) 同樣分布收斂於 F。

我們的目標是得到 T/√Var(T) = (Xbar - μ)/√Var(Xbar) 分布收斂於 N(0,1) 的結論，也就是有限群體之下中央極限定理成立。但 Xbar = Σ X_i 的成分 Xi 間並不獨立，即使把 Xbar 表示成 Σ Y_i t_i，諸 t_i 也相互不獨立，因此才費心構建另一樣本，另一隨機變數 T' = Σ (Yi - μ)I_[Ui≦f], 由於諸 Ui 相互獨立，諸指示函數 I_[Ui≦f] 也相互獨立。在 Yi 是系統性變化的情形，Yi 非隨機，而 T' 的諸項 (Yi - μ)I_[Ui≦f] 間也是相互獨立，因此要成立中央極限定理，我們需要的一個充分條件是 Linderberg 條件：

Σ ∫_[|Yi(I_[Ui≦f]-f)|≧ε√(nσ^2(1-f))] (Yi(I_[Ui≦f]-f))^2 dP/(nσ^2(1-f)) = 0

對任意 ε > 0 都成立。

我們也可直接看 T'/√Var(T') 的 m.g.f.

m(t) = Π_{i=1～N} [f e^{(1-f)Yi t/(√n σ)} + (1-f) e^{(-f)Yi t/(√n σ)}]

取對數（自然對數），則得

log(m(t)) = (t^2/2) f(1-f){ΣYi^2/[n(1-f)σ^2]}
+ (t^3/3!) f(1-f)(1-2f)ΣYi^3/[n(1-f)σ^2]^{3/2}
+ (t^4/4!) f(1-7f+12f-6f^3)ΣYi^4/[n(1-f)σ^2]^2 + ...
= t^2/2 + (t^3/3!){(1-2f)/√[n(1-f)]}(μ_3/σ^3)
+ (t^4/4!){(1-7f+12f-6f^3)/[n(1-f)^2]}(μ_4/σ^4) + ...

如果此處沒算錯，假設 0 << f << 1，假設在群體增大時 μ_k/σ^k 有界或增長的速度不如 n^{k/2-1}, k>2，則 n, N 同時趨於無窮而 f_N = n/N → f 時，T'/√Var(T') 的 m.g.f. 將以 e^{t^2/2} 為極限，也就是說 T'/√Var(T') 將分布收斂至標準常態，N(0,1)。在 N → ∞ 時，實務上假設 f < 1/2 是合理的，否則如果 f_N → 1，抽樣變成普查了；另一方面若 f_N → 0，實務上無妨假設先有 N → ∞ 然後才 n → ∞，則「有限群體」的問題消失，直接引用 i.i.d. 樣本的 CLT 即可。統計上引用 CLT 是實務問題，n 其實是有限，我們只要確定當前的群體、現在的抽樣，結果可以引用 CLT 做近似，因此本文假設 f_N = f 或 f_N → f 且 f 是小於 1 大於 0，而不考慮 n, N 的複雜行為。由前面 T 與 T' 的關係，我們得知 T/√Var(T) 有相同極限分布，亦即：有限群體抽樣之樣本平均數（或總值）亦適用中央極限定理。

關於有限群體抽樣之中央極限定理適用問題，前面的 Linderberg 條件看起來頗複雜，本文並未仔細討論；後面直接用動差母函數，但並未算出其對數展開式的一般項，而且要求群體動差系列的行為，也不很方便。不過，Hájek (1960) 就提出一個充分條件：

(max{(Yi - μ)^2}/σ^2)/min{n, N-n} → 0

Lehmann (1974) 的 Nonparametrics: Statistical Methods Based on Ranks 的充分條件是

[max{(Yi - μ)^2}/(Nσ^2)] max{(1-f)/f, f/(1-f)} → 0

其中 f 是 f_N = n/N 的簡寫，事實上此兩條件是等價的。較近的， P. K. Sen (1995) 也對 Hájek 的條件做了些討論。

我們以群體 Y = {1, 2, ..., N} 為例，μ = (N+1)/2, σ^2 = (N^2-1)/12，max{(Yi - μ)^2} = (N-1)^2/4，Hájek 條件是數學上考慮 f 因 N 之增大而可能趨近於 0 或 1，如果限制 f < 1/2, 也就相當於要求 f 趨近於 0 的速率必須慢於 O(1/N)，其實這也相當於 N → ∞ 時 n 也要趨於無窮。但 CLT 在數學上本來就要求 n → ∞，所以可以說在此群體架構之下，CLT 是可用的。基於 ranks 之非參數方法，假設原群體具連績型分布，對抽取之樣本觀測值都轉成 ranks，而後基於這些秩取某些部分做彙總，如兩獨立樣本之 Mann-Whitney 檢定，把兩樣本綜合排出 ranks，然後以一樣本之秩和為基礎建立檢定統計量，在虛無假說（兩群體相同分布）相當於從 N = n + m 個觀測值群體 {1, 2, ..., n+m} 中以不放還法隨機抽取一個大小為 n 的樣本（樣本 1 的大小是 n，樣本 2 是 m），依據有限群體之 CLT，若 N, n 夠大，標準化的秩和統計量接近常態分布。

在上面 ranks 群體的例子，σ^2 = (N^2-1)/12, 而群體之原點 k 階動差是 O(N^k)，所以中心動差也是 O(N^k) 或更低，前述動差序列條件成立，m(t) → e^{-t^2/2} 成立。另一個例子，假設 Y 群體是抽自一超群體，那麼 N 增大時只是從這個超群體隨機取得更多觀測值，因此群體 μ_k 在 N 增大時趨於穩定，所以動差系列的條件也滿足，CLT 再次成立。這也就是說：從有限群體抽出之樣本，不管是當做背後超群體的樣本，或看成是當前有限群體的樣本，都可適用 CLT，差別是：在觀念上是代表當前有限群體，或是代表背後的超群體？在數量上前者 Xbar 的變異數接近 (1-f)σ^2/n，後者則是 σ^2/n，因為考慮諸 Yi 的隨機性所以 Xbar 的不確定性提高。