樣本分位數的大樣本漸近分布－劉應興的部落格

假設群體分布是連續型，其支撑 (support) 是一個閉區間，在此區間有連續的 p.d.f. f(x)，故其分布函數 F(x) 在此區間是嚴格遞增的，對於任一 p, 0 < p < 1, 恰有一 x 值使 F(x) = p。令

g(p) = F^(-1)(p), 0 < p < 1

這樣的設定很嚴苛，只限於少數連續型機率分布，所幸我們常見的連續型分布就是這樣子的。

我們將樣本分位數定義為：將 X1, ..., Xn 由小而大排序成

X(1) ≦ . . . ≦ X(n)

後，取 k = {np}, 不小於 np 的最小整數，也可以表示成：

X({np}) = min {Xi, i = 1, ..., n : Σ_{j=1~n} I(Xj≦Xi) ≧ np}
= min {x in {X1, ..., Xn} : Σ_{j=1~n} I(Xj≦x) ≧ np}
= min {x : Σ_{j=1~n} I(Xj≦x) ≧ np}

式中 I(．) 是指標函數，符合括弧內條件則取值 1 否則 0。上列樣本 p 分位數的定義相當於經驗分布函數

Fn(x) = (1/n) Σ_j I(Xj ≦ x)

之虛擬反函數的定義

G(p) = inf{x; F(x) ≧ p} = sup{x; F(x) < p}, p in (0,1)

故亦可改寫樣本 p-分位數定義為

X({np}) = sup {x : Σ_{j=1~n} I(Xj≦x) < np}

但先前定義之另兩種寫法只在 Xi 間測試，就不適當了。

事實上樣本分位數可能有不同定義，例如取 X([np])（[np] 為不大於 np 的最大整數，或 (n+1)p 取代 np, 並且可能取 np 或 (n+1)p 最接近的兩樣本觀測值做簡單或加權平均。不過，不同定義之間不過是在 X([np]) 和

X({(n+1)p}) ≦ X({np}+1) ≦ X([np]+2)

之間取一值。當 n 很大時，其差距無限縮小。為什麼如此，因為 F 是連續型，F(X) 將在 (0,1) 區間形成均勻分布，也就是

F(X(1)) < F(X(2)) < . . . < F(X(n))

形成從 U(0,1) 抽取大小為 n 的樣本的順序統計量。因此，

我們先考慮 X 服從 U(0,1) 的情形。我們假設 j < k，令

U = X(k)-X(j), V = X(J) 其反變換 X(j) = V, X(k) = U+V

則 U, V 之聯合機率密度為

h(u,v) ={n!/[(j-1)!(k-j-1)!(n-k)!]}v^(j-1)u^(k-j-1)(1-u-v)^(n-k), 0 < v < u+v < 1

對 v 積分即是 u 的邊際 p.d.f.

g(u) = {n!/[(j-1)!(k-j-1)!(n-k)!]}u^(k-j-1)(1-u)^(n-k+j) ×
∫_[0,1-u] [v/(1-u)]^(j-1)[(1-u-v)/(1-u)]^(n-k)[1/(1-u)] dv
= {n!/[(j-1)!(k-j-1)!(n-k)!]}{(j-1)!(n-k)!/(n-k+j)!}u^(k-j-1)(1-u)^(n-k+j)
= {n!/[(k-j-1)!(n-k+j)!}u^(k-j-1)(1-u)^(n-k+j)

如果 k - j = r 為定數，則

P[U > ε] = ∫_[ε,1] {n!/[(r-1)!(n-r)!]} u^r (1-u)^(n-r) du
≦ ∫_[ε,1] {n!/[(r-1)!(n-r)!]} (1-u)^(n-r) du
= {n!/(r-1)!(n-r+1)!}(1-ε)^(n-r+1)

對任意 ε > 0，當 n → ∞ 時都得 P[U > ε] → 0，即 X(k)-X(k-r) 機率收斂於 0。不僅如此，

若 0 < α < 1，則 n^α U 機率收斂於 0。

因此至少我們能說：當群體分布是 U(0,1) 時，

X(k) = X(k-r) + o_p(1/√n)

也就是說 √n(X(k)-X(k-r)) → 0 in P, 當 n → ∞, r 固定。對一般連續型分布如本文開始的設定，則

F(X(1)) < F(X(2)) < . . . < F(X(n))

是取自 U(0,1) 之大小為 n 的隨機樣本順序統計量，因此

F(X(k)) = F(X(k-r)) + o_p(1/√n)

但我們假設 F 的支撑是一個區間，明確地說，是區間 (a,b) 的閉包 (closure), 即 (a,b) 本身加上端點，如果端點存在。Ｆ是嚴格遞增，連續，所以其反函數 g 也連續，

X(k) = g(F(X(k))) = g(F(X(k-r))) + o_p(1/√n) = X(k-r) + o_p(1/√n)

這也就是說，我們不管設定怎樣的規則從 X([np]) 到 X({(n+1)p}) 取樣本 p-分位數，並不影響其大樣本漸近分布（假設看的是它和某定點的離差放大為 √n 倍）。

令 k = {np}, 則

P[X(k) ≦ x] = P[k 個以上 Xi ≦ x] = P[Fn(x) ≧ k/n]

我們要的是 Zn = √n (X(k) - ξ) 的漸近分布。將 X(k) 表示成 ξ + Ｚn/√n。則

P[Zn ≦ z] = P[X(k) ≦ ξ + z/√n] = P[Fn(ξ + z/√n) ≧ k/n]

注意 np ≦ {np} = k ≦ np+1, 所以 √n(k/n - p) = O(1/√n) = o(1)。但

Fn(ξ + z/√n) = F(ξ + z/√n) + (1/n) Σ [I(Xi ≦ ξ + z/√n) - F(ξ + z/√n)]

令 Yni = I(Xi ≦ ξ + z/√n) - F(ξ + z/√n)，則

Yni = 1- F(ξ + z/√n) 機率 F(ξ + z/√n); = -F(ξ + z/√n) 機率 1 - F(ξ + z/√n)

故 E[Yni] = 0 而

E[Yni^2] = F(ξ + z/√n)(1- F(ξ + z/√n))
= (F(ξ)+f(ξ)z/√n+o(z/√n))[1-(F(ξ)+f(ξ)z/√n+o(z/√n))]
= p(1-p) + O(1/√n)

故

P[Fn(ξ + z/√n) ≧ k/n]
= P[Fn(ξ + z/√n) - F(ξ + z/√n) ≧ k/n - (F(ξ)+f(ξ)z/√n+o(z/√n)]
= P[(1/√n)ΣYni ≧ -f(ξ)z/√n + o(z/√n)]

然而 Yni, i = 1, ..., n 這個三角行列，每一橫列是 n 個 i.i.d. 隨機變數，均數 0 而變異數是 p(1-p)+O(1√n)，絕對符合 Lindeberg 條件，事實上它也滿足 Lyapounov 條件，或直接適用 i.i.d. 的情形，只是 Var(Yni) 多了個 O(1/√n)。總之，

(1/√n)ΣYni／√[p*(1-p*)] 向標準常態做分布收斂，

其中 p* = F(ξ + z/√n) = F(ξ) + f(ξ)z/√n + o(z/√n) = p + O(z/√n)　故

(1/√n)ΣYni／√[p(1-p)] 亦向標準常態做分布收斂。

所以 (1/√n)ΣYni 收斂於 N(0, p(1-p))，而

P[Fn(ξ + z/√n) ≧ k/n] → 1 - Φ(-f(ξ)z/√[p(1-p)]) = Φ(f(ξ)z/√[p(1-p)])

式中 Φ 為標準常態分布的分布函數。但上式左邊等於 P[√n (X(k) - ξ) ≦ z]，故

√n (X(k) - ξ) 向 N(0, p(1-p)/(f(ξ))^2) 做分布收斂

或者說：X({np}) 的漸近分布是 N(ξ, p(1-p)/[n(f(ξ))^2])。

如果同時考慮多個分位數，例如 0 < p < q < 1，令 Ynip 如前面 Yni 之定義，而 Yniq 是把 Ynip 的 p-分位數換成 q-分位數，則可證得 (1/√n) Σ_i (Ynip, Yniq) 收斂至常態，所以最後可推得 (X({np}), X({nq})) 的漸近分布也是常態，兩樣本分位數的漸近共變異數是 p(1-q)/[n f(ξp)f(ξq)]。擴充至多個分位數是很直接的，重點是樣本 p-, q-分位數的漸迎共變異數是

(min{p,q}-pq)/[n f(ξp)f(ξq)]

如果是多變量群體，例如雙變量，(Xi,Yi) 是隨機樣本，分位數不論群體或樣本都是各變數分別定義，則再一次可得 (X({np}),Y({nq})) 的聯合分布漸近多（雙）變量常態，其漸近共變異數是

(F(ξp1,,ξq2) - pq)/f(ξp1,,ξq2)

多變量群體和雙變量群體沒有根本的不同。

如果群體分布不像設定的那樣，例如 F(x) = p 不是唯一解，分位數雖仍可唯一定義，卻不是根本上的唯一。並且，此時意謂在某一區間 f(x) = 0，而在所定義的群體分位數

ξp = inf{x; F(x) ≧ p}

這一點 f 的定義其實是隨意的，它可能是 F 在此點的左導數，而右導數是 0。如果 f(ξ) = 0, 上面的「漸近分布」就無意義了，因為漸近變異數分母是 0。如果 F 在 ξ 不可微呢？上面有些時候要做 F 在 ξ 的線性近似，不可微是不行的。如果雙邊的單邊導數都存在但不等，則因

P[X(k) ≦ ξ + z/√n] = P[Fn(ξ + z/√n) ≧ k/n] → Φ(f(ξ)z/√[p(1-p)])

當 z > 0 時 F 在右邊對對 ξ 展開；z < 0 時 F 在左邊展開，因此 X({np}) 的漸近分布將是兩個半常態分布拼湊而成，在 ξ 右邊用 F 的右邊導數；在 ξ 左邊用 F 的左邊導數。

有另一想法考慮

Yn = Σ_{i = 1~n} I(Xi ≦ ξ) ～ binom(n,F(ξ))

換個方式表示，Yn = n Fn(ξ)，其大樣本漸近分布為 N(n F(ξ), n F(ξ)(1-F(ξ))，試圖由此來導出樣本 p-分位數，例如 X({np}) 的漸近分布。

取反函數變換 g(u) = F^(-1)(u), 則依 delta-method，

g(Yn/n) asymptotic distributed as N(ξ, (g'(F(ξ)))^2F(ξ)(1-F(ξ))/n)

也就是

√n(g(Yn/n) - ξ) converges in distribution to N(0, F(ξ)(1-F(ξ)/(n (f(ξ))^2))

其中 g'(F(ξ)) = 1/F'(g(F(ξ))) = 1/f(ξ)。所謂 delta-method 就是將 g(Yn/n) 在 p 點展開，

g(Yn/n) - g(p) = g'(p)(Y/n - p) + o(Yn/n-p) = g'(p)(Y/n - p) + o_p(1/√n)

所以 Yn 的漸近分布被轉成 g(Yn/n) 的漸近分布。但這和我們所要的 X({np}) 的漸近分布有什麼關係？

首先我們來看看 X({np}), ξ, p 和 Yn/n 的關係：假設 Fn 定義了前述的虛擬反函數 g*，使得

Fn(g*(Fn(x)) = x for all real x

則

g*(p) = X({np}), Yn/n = Fn(ξ), g(p) = ξ, p = F(ξ)

樣本分位數的大樣本漸近分布

因為 F 是連續型，假設樣本觀測值互不相等，即

X(1) < X(2) < ... < X(n)

所以 X({np}) 與 ξ 之間所包含的樣本觀測值個數，等同於 Fn(Xi) 在 p 與 Yn/n 之間所包含的觀測值個數，因為它們就是在 Fn 曲線的 A, B 兩點之間所含的觀測值個數。因為假設 Xi 間互不相等，|Yn-np| 與此個數之差頂多是 1，

Fn(X({np})) - p = {np}/n - p = - (Yn/n - p) + O(1/n)

其中 O(1/n) 表示可以用 M/n 界定其大小，實際上根據上面結果，這裡 M 可取 1。

但是 A, B 兩點從橫向連接 F 曲線是高度 p 與 Yn/n，從縱向連接 F 曲線是卻是橫軸上的 X({np}) 與 ξ。我們要的是後者，也就是

X({np}) 或 X({np}) - ξ 的漸近分布

或其對應 F 的高度 F(X({np})) 或高度差 F(X({np})) - p。問題是 F(X({np}) 與 p 的距離是否也以 o_p(1/√n) 近似等於 Yn/n 與 p 的距離？或說：是否

F(X({np})) - F(ξ) = -(Fn(ξ) - Fn(X({np}))) +o_p(1/√n) ?

或者，等價地：能否證明

(Fn(X({np})) - F(X({np}))) - (Fn(ξ) - F(ξ)) = o_p(1/√n)

我們知道在每一點 x，

√n(Fn(x) - F(x)) → N(0. F(x)(1-F(x))) in d.（分布收斂）

或 √n(Fn(x) - F(x))／√[F(x)(1-F(x))] 向標準常態做分布收斂，因此，可以表示成

Fn(x) = F(x) + z(x)√{[F(x)(1-F(x))]/n} + o_p(1/√n)

將 X({np}) 寫成 X(k)，則

(Fn(X(k)) - F(X(k))) - (Fn(ξ) - F(ξ))
　　 = z(X(k))√{[F(X(k))(1-F(X(k)))]/n} - z(ξ)√{p(1-p)/n} + o_p(1/√n)

可證明 X(k) → ξ, z(x) 只是 Fn(x) - F(x) 的標準化表示，因此，上列右式前兩項的差應該可以用 o_p(1/√n) 表示，則

F(X({np})) - F(ξ) = -(Fn(ξ) - Fn(X({np}))) +o_p(1/√n) = -(Yn/n - p) + o_p(1/√n)

成立。所以 √n(F(X({np})) -p) 向 N(0, p(1-p)) 做分布收斂。同時，

√n (X(k) - ξ) = √n (g(F(X(k))) - g(F(ξ))) = -√n(g(Yn/n) - g(p)) + o_p(1)

所以 √n (X(k) - ξ) = -g'(ξ) √n (Yn/n - p) + o_p(1)。

也有作者討論離散型分布群體抽樣之順序統計量的問題，例如由 Yanyuan Ma, Marc G. Genton and Emanuel Parzen (2011) 所發表的「中點分布函數法」，將分布函數調整為

F*(x) = F(x) - (1/2)P[X = x] 當 P[X=x] > 0

在 F 的跳躍點，F(x) 由本來的右極限 F(x+) 調整為 (F(x-)+F(x+))/2, 左極限和右極限的中點，然後用直線段連接 F 各跳躍點調整過的 F(x) 值。如此一來，F*(x) 成為嚴格遞增函數；除最小點和最大點以外，處處存在左右導數。事實上如二項、Poisson 等實務常見的離散群體，可以採用下列方法做機率值連續化：

P[X = k] = p_k → 改為 f(x) = p_k, k-1/2 < x < k+1/2

此法也使得分布函數在一定範圍是連績的、嚴格遞增的、處處存在左右導數。不過 mid-distribution function 的方法左右導數不等的點在原 F 的跳躍點，而機率連續化法則是在各連續整數點的中點。似此類方法至少有兩個問題：一是改變了群體分位數的定義，且樣本觀測值僅能取自原分布之跳躍點；其二是在左右導數不等的點所對應的分位數，漸近分布並非常態，而是兩半常態曲線拼湊的分布（如前述）。事實上，離散群體的推論問題，是否需要中位數及分位數的推論及其大樣本性質，或有待商榷。

在絕對離差和準則下 Σ|Xi - a|，資料 X1, ..., Xn 的位置代表值 a 取樣本中位數 Md 最理想；對正負離差取不同權量，則 a 的最適值成為分位數，如「分位數迴歸」一文所述。雖然樣本分位數在某種準則下最適合資料，但想用它可用來估計群體分位數，卻只算是一種合理的想法，是資料配適 (data fitting) 的方法；而不是基於某種「估計量評估標準」，例如某種損失函數下期望損失，頻率論者的風險函數法或貝氏學派的後驗期望損失最小。也就是說：資料配適好的估計並不一定是一般所期望的最佳估計。像最小平方估計在標準線性模型、在平方誤差損失之下是線性估計中最佳，只是個巧合。怎麼這麼說？假設要估計 θ，採用某種資料配適準則

optimize K(θ; x) ==> 得 θ 之估計 T(x) = arg opt K(θ; x)

而統計決策理論的關點，這裡採頻率論觀點，是要

minimize R(δ, θ) = E[L(δ(x),θ)] for some meanning

假設 K(θ; x) 也是要找極小，它代表的是諸資料點 Xi 與 T(x) 的一種總和距離。把 T(x), θ 和資料點 x 都放在同一個空間，R^n，　則 T(x), θ 在 R^n 的「均數軸」X1 = ... = Xn 上。對每一個 x 我們在均數軸上找到最「接近」的點 T(x)；至於它和 θ 是否接近，根本管不著。在平方離差和準則下，T(x) = Xbar，而

K(θ; x) = ||θJ - x||^2 = n(θ-T(x))^2 + ||T(x)J - x||^2

第一項在線性、不偏的要求下期望值極小化 E[(θ-T(x))^2; θ] 完全可說是巧合。如果換成絕對離差和準則加絕對誤差損失，則

K(θ; x) = |θJ - x| ≧ n|θ-T(x)| + |T(x)J - x|

其中 |．| 在向量代表 1-norm（絕對值之和），在實數則是普通絕對值。此時有幾個問題：一是兩邊是不等式——這可能不重要，只是我們未能建立一個等式罷了，但它並不能說明什麼；其次，T(x) 取資料 x 的中位數使 |T(x)J - x| 最小，但對 |θ-T(x)| 似乎沒有直接的辦法可以說樣本中位數是個優良選擇，即使 θ 是群體中位數。

樣本中位數做為群體中位數的估計，究竟效率如何？比之其他估計量如樣本平均數呢？要比較必須是兩者估計的東西相同，所以我們限制在對稱分布群體。再者，有限樣本方法，或稱小樣本方法，或稱正確法 (exact approach)，是比較麻煩的。因此，這裡考慮大樣本方法，或稱漸近方法 (asymptotic approach)。首先是常態群體，樣本中位數漸近服從常態而樣本平均數一直是常態：

樣本中位數～. N(μ, (π/2)σ^2/n); 樣本平均數～ N(μ, σ^2/n)

不管比較兩者抽樣分布的標準差或平均差 E[|T(x)-μ|]，其漸近效率都是 √(π/2) = 0.798。Cauchy 分布 n 個觀測值的樣本平均數抽樣分布同於 n = 1 的分布，而且期望值不存在，平均誤差也是無窮大；但樣本中位數仍漸近服從 N(μ, (πσ)^2/(4n))。邏輯斯 (logistic) 分布可說相當接近常態分布，雖然函數形式不同，群體變異數 (πσ)^2/3，在中心點 μ 之密度為 1/(4σ)，所以大樣本

樣本中位數～. N(μ, 4σ^2/n); 樣本平均數～. N(μ, π^2σ^2/(3n))

樣本中位數漸近效率是 0.822。拉不拉斯 (Laplace) 分布，p.d.f. 是

f(x) = e^{-|x-μ|/σ}/(2σ), -∞ < x < ∞

樣本數夠大時得

樣本中位數～. N(μ, σ^2/n); 樣本平均數～. N(μ, 2σ^2/n)

立得樣本中位數漸近效率為 √2 = 1.414。均勻分布 f(x) = 1/σ, μ-σ/2 < x < μ+σ/2,

樣本中位數～. N(μ, σ^2/(4n)); 樣本平均數～. N(μ, σ^2/(12n))

樣本中位數相對於樣本平均數效率僅 1/√3 = 0.577。

以上是大樣本方法。在大樣本方法之下，幾種分布由於漸近分布是常態（除 Cauchy 分布群體之樣本平均數完全不收斂以外），因此絕對誤差損失和平方誤差損失的相對效率是相同的。至於小樣本方法，Ryszard Zieli´nski (2004) "Optimal Quantile Estimators Small Sample Approach" 考慮 L-估計量，即：以順序統計量之線性組合構建的估計量，討論了其中的一致最小變異不偏估計、一致最小絕對離差估計及其他評估標準下的最適估計等。

｛附｝

樣本分位數機率收斂至群體分位數

令 k = {np}, 前文有

P[X(k) ≦ x] = P[k 個以上 Xi ≦ x] = P[Fn(x) ≧ k/n]

取 x = ξ+ ε，則

P[X(k) ≦ ξ + ε] = P[Fn(ξ+ε) ≧ k/n] = P[Fn(ξ+ε) ≧ {np}/n]

依大數法則，Fn(x) → F(x)，再令 p(ε) = F(ξ+ε)，則

P[X(k) ≦ ξ + ε] = P[Fn(ξ+ε) - p(ε) ≧ k/n - p(ε)]

式中 Fn(ξ+ε) - p(ε) 當 n → ∞ 時向 0 做機率 1 收斂，當然也做機率收斂，故

lim P[X(k) ≦ ξ + ε] = 1 若 ε > 0; = 0 若 ε < 0

更嚴謹地說，對任意 δ > 0, γ > 0, 存在 N, 使得：

當 n > N 則 P[|Fn(ξ+ε) - p(ε)| > γ] < δ