原子核物理学においては固有値の並びの順序を考慮に入れた統計量が必要になるのですが,今回のコラムでは順序統計量の分布についての統計学的な性質をおさえておくことにしましょう.
累積分布関数をF(x),確率密度関数をf(x)とするある連続分布からのn個の観測値x1〜xnを大きさの順に並び換えたものを
x(1),・・・,x(n)
と書くことにします.x(1)は最小値,x(n)は最大値であり,x(k)をk番目の順序統計量といいます.
===================================
【1】最大値,最小値の分布
最大値x(n)の累積分布関数G(x)は
x(i)≦x (i=1〜n)
となること,すなわち,x1,・・・,xnのすべてがx以下になることと同値ですから,
P{x(1)<=x and x(2)<=x and x(n)<=x}
=P{x(1)<=x}P{x(2)<=x}・・・P{x(n)<=x}
=F(x)*F(x)*・・・*F(x)
より,
G(x)={F(x)}^n
したがって,最大値の確率密度関数g(x)はxについて微分して
g(x)=n{F(x)}^(n-1)f(x)
同様にして最小値x(1)の累積分布関数は
G(x)=1-{1-F(x)}^n
確率密度関数は
g(x)=n{1-F(x)}^(n-1)f(x)
となります.
一般に,順序統計量x(k)の確率密度関数は,x(k)≦xとなることとx1,・・・,xnのうちk個以上がx以下になることと同値であり,
P{x(k)<=x}=ΣnCkF(x)^k{1-F(x)}^(n-k)
より
g(x)=n!/(k-1)!(n-k)!F(x)^(k-1){1-F(x)}^(n-k)f(x)
で与えられます.
n!/(k-1)!(n-k)!
は2項係数nCk=n!/k!(n-k)!を用いてknCk,ベータ関数を用いて1/B(k,n-k+1)と表せますから
g(x)=knCkF(x)^(k-1){1-F(x)}^(n-k)f(x)
g(x)=1/B(k,n-k+1)F(x)^(k-1){1-F(x)}^(n-k)f(x)
などとも書くことができます.
[1]x1,・・・,xnが互いに独立に区間(0,1)の一様分布:F(x)=xに従っているとき,x(k)の確率密度関数は
g(x)=1/B(k,n-k+1)x^(k-1){1-x}^(n-k)
すなわちベータ分布:Beta(k,n-k+1)であることがわかります.
したがって,
E[x(k)]=k/(n+1)
V[x(k)]=k(n-k+1)/(n+1)^2/(n+2)
となります.最大値の分布に関して
E[x(n)]=n/(n+1)
ですから,nが無限に大きくなるとx(n)は1に確率収束することを意味しています.これは直感的にも自明です.
[2]次に,2変数の場合で考えて見ましょう.x,yがパラメータλの指数分布
F(x)=1−exp(-λx)
にしたがい,互いに独立なとき,最小値の累積分布関数は
1-{exp(-λx)}^2 =1-exp(-2λx)
ですから,最小値の分布はパラメータ2λの指数分布になります.
n変数の場合の最小値の分布関数は
1-exp(-nλx)
すなわちパラメータnλの指数分布になります(平均:1/nλ,分散:1/(nλ)^2).しかし,最大値の分布は{1-exp(-λx)}^nですから指数分布にはなりません.
具体的な分布型は煩雑になりますが,指数分布の順序統計量の平均・分散に関しては,それぞれ引数が整数のときのジガンマ関数・トリガンマ関数に帰着し,簡単な形になります.
E[x(k)]=1/λΣ1/(n-k+1)
V[x(k)]=1/λ^2Σ1/(n-k+1)^2
したがって,n→∞のとき,最大値x(n)の分布の平均は無限に大きくなりますが,ほぼ
1/λ・(logn+γ) (γはオイラーの定数:0.57722・・・)
に等くなります.また,Σ1/k^2=π^2/6=ζ(2)ですから,分散は
1/λ^2・π^2/6
に収束します.
[3]正規分布:
標準正規分布
f(x)=1/(2π)^(1/2)exp(-x^2/2)
の場合,
n=2:E[x(2)]=-E[x(1)]=1/π^(1/2)
n=3:E[x(3)]=-E[x(1)]=3/2π^(1/2),E[x(2)}=0
n=4:E[x(4)]=-E[x(1)]=6arctan(2^(1/2))/π^(3/2)
E[x(3)}=-E[x(2)]=6/π^(1/2)-18arctan(2^(1/2))/π^(3/2)0
n=5:E[x(5)]=-E[x(1)]=15arctan(2^(1/2))/π^(3/2)-5/2π^(1/2)
E[x(4)}=-E[x(2)]=10/π^(1/2)-30arctan(2^(1/2))/π^(3/2)
E[x(3)]=0
一般に
φ(x)=1/(2π)^(1/2)σexp(-(x-μ)^2/2σ^2)
の場合,最大値の分布は
G(x)={Φ((x-μ)/σ)}^n
g(x)=n{Φ((x-μ)/σ)}^(n-1)φ((x-μ)/σ)/σ
E[x(n)]=μ+e1σ (e1はnについての増加関数)
V[x(n)]=e2σ^2 (e2はnについての減少関数)
最小値の分布は
G(x)=1-{1-Φ((x-μ)/σ)}^n
g(x)=n{1-Φ((x-μ)/σ)}^(n-1)φ((x-μ)/σ)/σ
E[x(1)]=μ-e1σ (e1はnについての増加関数)
V[x(1)]=e2σ^2 (e2はnについての減少関数)
と表すことができます.
[4]半円分布:f(x)=(1-x^2)^(1/2) (-1≦x≦1)の場合,
F(x)=1/2[x(1-x^2)+arcsinx]+π/4
ですから,最大値x(n)の累積分布関数G(x)は
G(x)={1/2[x(1-x^2)+arcsinx]+π/4}^n
となり,具体的な分布型は煩雑です.これについては(その2)で述べることにします.
===================================
【2】極値の極限分布
最大値・最小値などの極値の分布は,有限のnについては【1】で与えられましたが,これらの分布はnとともにずれていきます.
そこで,一様分布の最大値[1]において,
y(n)=(n+1){x(n)-n/(n+1)}
と変数変換すると,その分布関数はF(x(n))=x^nより
G(y(n))={(y+n)/(n+1)}^n={1+(y-1)/(n+1)}^n
ですから,n→∞のとき,その極限分布は指数分布
G(y)=exp(-(1-y))
になることが理解されます.
また,指数分布の最大値[2]では,
y(n)=x(n)-logn
とおくと,
G(y(n))={F(y+logn)}^n={1-exp(-y)/n}^n
ですから,n→∞のとき,その極限分布は
G(y)=exp(-exp(-y))
すなわち2重指数分布の形に表されます.
E[y]=γ (γはオイラーの定数:0.57722・・・)
V[y]=π^2/6=ζ(2)
このようにnが大きいとき,漸近的にある位置=尺度モデルに収束するための条件
{G(x)}^n=G(a+bx)
の解を考えると,n→∞としたときの極値極限分布は,ワイブル分布(指数分布を含む)か2重指数分布(ガンベル分布)の2つのうちどちらかに限られることが証明されています.ワイブル分布は指数分布にしたがう確率変数のベキ乗変換であり,一方,2重指数分布は指数分布にしたがう確率変数の対数変換として導かれる分布です.
一般の分布について,その極限分布が上にあげた2つの型のいずれかになるというのは非常に有益な結果であり,しかも独立同分布の仮定が多少くずれたとしてもワイブル分布や2重指数分布があてはまることが多いとされています.
また,分布の裾が軽い場合の最大値の極限分布としては2重指数分布が現れることも知られていて,たとえば,正規分布[3]
f(x)=1/(2π)^(1/2)exp(-x^2/2)
については,
y(n)=(2logn)^(1/2){x(n)-(2logn)^(1/2)}
の分布が2重指数分布に近づくことが示されます.
===================================
【3】範囲の分布
母集団からn個の標本を抽出したとき,最大値x(n)と最小値x(1)の差
r=x(n)-x(1)
を範囲といってばらつきの程度を示す1つの指標となります.ここでは範囲の確率分布を示しておきます.
G(r)=n∫(-∞〜∞)[F(x+r)-F(x)]^(n-1)f(x)dx
g(r)=n(n-1)∫(-∞〜∞)[F(x+r)-F(x)]^(n-2)f(x)f(x+r)dx
[1]一様分布:f(x)=1/θ 0≦x≦θ F(x)=x/θ
g(r)=n(n-1)∫(0,θ)(r/θ)^(n-2)(1/θ)^2dx=n(n-1)r^(n-2)/θ^n
[2]指数分布:f(x)=θexp(-θx) F(x)=1-exp(-θx)
g(r)=n(n-1)∫(0,∞)(exp(-θx)-exp(-θ(x+r)))^(n-2)θexp(-θx)θexp(-θ(x+r))dx=(n-1)exp(-θr){1-exp(-θr)}^(n-2)/θ
母集団分布が一様分布や指数分布のとき,これは簡単な形になりますが,正規母集団のときには簡単にはなりません.正規分布の場合の期待値および分散は
E[r]=d1σ (d1=2e1)
V[r]=d2σ^2
なお,2重指数分布(ガンベル分布)の範囲の分布については第2種変形ベッセル関数を用いて表されることを申し添えておきます.ミッドレンジ:{x(1)+x(n)}/2の分布については省略します.
===================================
【4】標本中央値の分布
簡単のために,標本の大きさnを奇数とし,n=2m+1とおきます.【1】において,n=2m+1のとき,中央値x(m+1)の確率密度関数は
g(x)=(2m+1)!/(m!)^2F(x)^m{1-F(x)}^mf(x)
で与えられます.母集団分布が一様分布や指数分布のときこれは簡単な形になりますが,正規母集団のときには簡単にはなりません.
[1]確率変数が区間(0,1)の一様分布にしたがうとき,標本中央値x(m+1)については,
E[x(m+1)]=1/2
V[x(m+1)]=1/4(2m+3)=1/4(n+2)
一様分布の標本平均値と平均値の分散はそれぞれ1/2,1/(12n)ですから,標本平均値と標本中央値は一致しますが,その分散はn>1のときには常に標本中央値のほうが大きくなります.
[2]コーシー分布の標本中央値
コーシー分布
f(x)=1/π・α/{α^2+(x-μ)^2}
F(x)=1/π[arctan(x-μ)/α]+1/2
より,その分布は
g(x)=(2m+1)!/(m!)^2π2^2m{1-(2/πarctanx-μ)/α)^2}^mα/{α^2+(x-μ)^2}
となります.
長い積分計算の後,
期待値:E[x(m+1)]=μ
分散 :V[x(m+1)]=α^2/(n+2)(π/2)^2{1+2/(n+4)(π/2)^2+3/(n+4)(n+6)(π/2)^4+・・・}
したがって,nが十分大きいところでは
V[x(m+1)]=α^2/(n+2)(π/2)^2
これにより,標本中央値の分散は標本の大きさnを大きくすると小さくなることが示されました.コーシー分布に従う変数については,標本平均値に関する中心極限定理が成り立たないわけですから,まことに注目すべきことです.
===================================
【5】標本中央値の漸近分布
正規母集団以外であっても,その母平均,母分散をそれぞれμ,σ^2として,標本平均値の分布が漸近的に正規分布N(μ,σ^2/n)になることは高校の教科書にも取り上げられて,標本平均値についての統計的性質は「中心極限定理」としてよく知られています.しかし,標本中央値の漸近分布を取り上げたものは少ないようです.
以下では,標本中央値に関する極限定理「母集団のメジアンをμmとすると,メジアンの分布は漸近的に正規分布N(μm,1/{4n[f(μm)]^2})になる」ことを証明していきます.
(証明)xがg(x)=(2m+1)!/(m!)^2F(x)^m{1-F(x)}^mf(x)にしたがうとき,
u=(x-μm)/(1/{4n[f(μm)]^2})^(1/2)=2{nf(μm)(x-μm)}^(1/2)
が漸近的にN(0,1)にしたがうことを示せばよい.
uの確率密度関数は,
x=u/2{nf(μm)+μm}^(1/2),dx=du/2{nf(μm)}^(1/2)
より,
h(u)du=g(x)dx=g(u/2{nf(μm)+μm){^(1/2)/2{nf(μm)}^(1/2)du
h(u)=g(u/2{nf(μm)+μm)}^(1/2)/2{nf(μm)}^(1/2)
=(2m+1)!/(m!)2*{∫(-∞,u)f(t)dt*∫(u,∞)f(t)dt}^mf(u/2{nf(μm)+μm)}^(1/2)/2{nf(μm)}^(1/2)
ここで,n→∞のときの極限を考える.スターリングの法則
m!=(2π)^(1/2)m^(m+1/2)exp(-m)
(2m+1)!=(2π)^(1/2)(2m+1)^(2m+3/2)exp(-2m-1)
を使って簡約化すると
(2m+1)!/(m!)2→1/(2π)^(1/2)*{n*2^n}^(1/2)
また,積分学における第1平均値定理により
{∫(-∞,u)f(t)dt*∫(u,∞)f(t)dt}^m
=2^(-2m){1-[uf(ξ)]^2/n[f(μm)]^2}^m → 2/2^nexp(-u^2/2)
f(u/2{nf(μm)+μm)}^(1/2)/2{nf(μm)}^(1/2) → 1/2(n)^(1/2)
したがって,
h(u)→1/(2π)^(1/2)exp(-u^2/2) 〜 N(0,1)
[補]スターリングの公式
n!=(2π)^(1/2)n^(n+1/2)exp(-n)
より
2nCn 〜 2^(2n)/(πn)^(1/2)
−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
さらに,F(q)=pなる統計量qの漸近分布は
N(q,p(1-p)/{n[f(μm)]^2})
となることを示すこともできます.この式でp=1/2とおくと中央値の漸近分布,p=1/4とおくと第1四分位数の漸近分布が得られます.また,順序統計量の線形結合の分布に対しても,漸近正規性が成り立つことが証明されています.
母集団分布が正規分布N(μ,σ^2)のとき,
f(x)=1/(2π)^(1/2)σexp{-(x-μ)^2/2σ^2}
ですから,
μm=μ
f(μm)=1/(2π)^(1/2)σ
となるので,標本中央値の分布は漸近的にN(μ,πσ^2/2n)になります.すなわち,母集団分布が正規分布に従うとき,標本平均値と標本中央値は一致しますが,分散は標本中央値のほうが大きくなり,標本中央値の分散は標本平均値の分散のπ/2倍になることが示されます.
これより,標本中央値は標本平均に比べて
(σ^2/n)/(πσ^2/2n)=2/π=0.637(約63.7%)
の有効性しかもたないことがわかります.これを漸近相対効率(ARE:asymptotic relative efficiency)といいます.
また,コーシー分布
f(x)=1/π・α/{α^2+(x-μ)^2}
の中央値の分布が漸近的に平均μ,分散π^2α^2/4nの正規分布になることも導き出されます.これは【4】の結果:α^2/(n+2)(π/2)^2と漸近的に等しくなります.したがって,
N(μm,1/{4n[f(μm)]^2})
式は簡単な推定方式ながら,かなりよい推定量を与えてくれることがわかります.(必ずしも極限定理の分散のほうが小さいというわけではありません.)
===================================
【6】標本中央値のロバスト性
標本平均は,母集団が正規分布N(μ,σ^2)に従うとき,それぞれ母平均μの不偏推定量でありかつあらゆる不偏推定量のなかで分散が最小です.したがって,正規分布の平均を推定するには標本平均が一番よい推定量です.しかしながら,母分布が正規でなかったりという状況では必ずしも最適な推定量とはいえませんし,また,標本平均は外れ値の影響を受けやすいという欠点もあります.
標本平均のこの欠点を標本中央値に対する漸近相対効率(ARE)を調べることによって見てみます.前節で示したように
ARE=(中央値の漸近分散/標本平均の漸近分散)^(-1)=4f^2(0)σ^2
で定義されます.
たとえば,自由度νのt分布であれば母分散はν/(ν-2)ですから,標本平均の漸近分散はν/(ν-2)/n,一方,標本中央値の漸近分散は
1/4n{f(0)}^2=1/4νπΓ^2(ν/2)/Γ^2((ν+1)/2)/4n
ですから,まとめると以下の表が得られます.
ARE
正規分布 0.64
t分布 (df=3) 1.62
(df=4) 1.13
(df=5) 0.96
(df=10) 0.76
両側指数分布 2
ロジスティック分布 0.82
コーシー分布 ∞
正規分布の場合は標本中央値の分散のほうが大きくなりましたが,表からわかるように,正規分布より少しでも裾が長いと想像されるときには,それがどんな分布であっても,位置母数推定に関しては,標本中央値のほうが標本平均よりAREが1より大きいか,または小さくてもそれ程小さくないという点で優れていると考えられます.この性質を称して,標本中央値は位置母数のロバスト推定量であるといいます.
一般的にいって,順序統計量は外れ値や分布形の非対称性の影響を受けにくい指標であり,たとえば,平均値と中央値を比べてみると他のデータと非常にかけ離れた値が混入している場合でも中央値のほうが平均値より狭い範囲に密集するため,安定した結果を与えてくれます.このように極端な値が混入しても推定量が大きく変動することのないような性質を持った推定量をロバスト推定量といいます.
標本平均や標本分散は外れ値の影響を受けやすく頑健とはいえませんが,標本中央値や標本4分位数は外れ値に対して抵抗性があります.そのため,中心位置のロバスト推定量としては中央値,散らばりのロバスト推定量としては4分位偏差や範囲がよく用いられます.正規分布より相当にばらつきに大きいことが想定される場合かあるいは分布形に関する情報が乏しいときは,ロバスト推定量の使用が薦められます.
===================================