■どの確率モデルを選定するか(その4)

 正にゆがんだ分布でも負にゆがんだ分布でも,適用できる分布としてベータ分布があげられますが,使い勝手のいい分布として,ラムダ分布なるものがあるそうです.

 その累積分位関数と分位密度関数は,位置母数をλ1,尺度母数をλ2,形状母数をλ3,λ4として,

  F(x)=λ1+{x^λ3−(1−x)^λ4}/λ2

  f(x)={λ3x^λ3-1+λ4(1−x)^λ4-1}/λ2

  0<x<1

で与えられます.

 かなり前のことですが,ヒストグラムに対する分布関数の適合プログラムを作成したことがありました.連続分布,離散分分布併せて100種類もの分布関数を扱うものでしたが,その中ではサポートされておらず,私にとっては初めて知ることになった分布です.

 その詳細は

  [参]四辻哲章「数値データ適合分布」プレアデス出版

に譲ることにして,今回のコラムでは順序統計量について考えてみることにします.

===================================

[1]ロバスト推定法

 誤差の分布が正規分布の場合だけではなく,裾の広い非正規分布に関してもパラメータ推定値の誤差分散が十分小さくなるような推定法が望ましいわけで,このような推定法をロバスト推定法といいます.ロバスト推定には大別して,M推定(最尤法に基づくロバスト推定),L推定(順序統計量に基づくロバスト推定),R推定(順位統計量に基づくロバスト推定)がありますが,M推定法の名前は最尤法(maximum likelihood method)から拡張されたということに由来しています.Lは線形結合(linear combination),Rは順位(rank statistic)の頭文字にちなんでいます.

===================================

[2]順序統計量と順位統計量

 母分布が正規分布のときの検定や推定がパラメトリック論ですが,データノ中に異常値(外れ値)を含む場合はデータ解析の結果が異常値に振り回されやすい,裾の重い分布のときの推測法しとて良くないなどの欠点があります.

 約100年に及ぶ近代統計学の歴史の中で,標本の母集団として正規分布を仮定した統計学(正規分布統計学)の時代を経て,順序統計学や順位統計学などいろいろな考えに基づく種々の統計学が誕生しました.これらを別々の統計学としてそのままにせず,ロバストネスの観点で統合することを試みてみます.

===================================

[3]標本分布の五数要約(five number summary)

 観測データの解析は度数分布図を描くことから始まります.記述統計学においてはその分布の形についての情報をなるべく失わないようにデータを要約することが重要です.

 データを統計パッケージに入力すると,平均値,分散,標準偏差,変動係数,歪度,尖度・・・などいろいろな統計量が出力されてきますが,この中で,データの位置情報に対応する指標としての平均値と,尺度情報(散らばり)に対応する指標としての標準偏差・分散がデータ全体を要約する統計量として検定・推定統計学にとって最も重要な統計量であることはいうまでもありません.

 そのため,データの分布形はしばしばただ2つの統計量,「平均値±標準偏差」で要約され,正規分布であればこの範囲内に入るデータの割合は約68.3%となり68%信頼区間を表わすと見なすことができます.しかし,平均値・分散・標準偏差だけでは左右対称でない分布や外れ値を含む場合に要約統計量とすることには問題があります.実は,データを平均値と標準偏差で要約しようとするときには暗黙のうちに正規分布もしくはそれに近い分布を想定しているのです.

 データを整理して図表にまとめる場合に多くの論文で「平均値±標準偏差」に要約して記載されていますが,要約の過程で分布に関する情報の一部が切り捨てられるので,裾の長い非対称な分布形の持つデータでは「平均値±標準偏差」で表現することは不適切であるばかりか,正規性・等分散性の仮定を無視することからパラメトリック統計処理で誤った結論を導く可能性が大きくなります.

 そこで,データ分布が非対称形を示す場合においてもデータをより適切に要約する非常に簡単な方法として,データを大きさの順に並べ換えたもの(順序統計量)から0%,25%,50%,75%,100%の位置のデータを要約値として採用する方法が考え出されています.これを五数要約といい,それぞれ最小値,第1四分位数,中央値(メディアン),第3四分位数,最大値に対応します.五数要約では,「平均値±標準偏差」の代わりに,代表的な三数を用いて,しばしば50%[25%,75%]と表現され,[ ]の範囲内にデータの50%が存在するという分布状況が正確に把握できます.

 四分位統計量を使った五数要約データをグラフ表示し,その分布傾向をよりいっそう明確かつ直感的に把握することができるように工夫した作図法に箱ひげ図があり,簡単なグラフのわりには多くの情報を含んでいるため,最近よく使われています.この表記法の利点は位置情報,ばらつき情報のみならず分布の形に関しての情報も得ることができます.とくに,裾が長い分布をもつデータにおいてはデータの分布状況が正確に把握できることにあります.また,範囲(最大値と最小値の差)や4分位偏差(第3四分位数と第1四分位数の差)なども散らばり方の目安としてよく用いられます.

===================================

[4]順序統計量の分布

 五数要約では順序統計量を扱いましたが,順序統計量の分布についておさえておくことにしましょう.累積分布関数を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-1f(x)

 同様にして最小値x(1)の累積分布関数は

  G(x)=1-{1-F(x)}^n

確率密度関数は

  g(x)=n{1-F(x)}^n-1f(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)

すなわちベータ分布であることがわかります.したがって,

  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+γ)(γはオイラーの定数:.577・・・)に等くなります.また,(Σ1/k2=π2/6=ζ(2))ですから,分散は1/λ^2π^2/6に収束します.

−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

(2)極値の極限分布

 最大値・最小値などの極値の分布は有限のnについては前述の式で与えられましたが,これらの分布は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重指数分布の形に表されます.

 このように,nが大きいとき,漸近的にある位置=尺度モデルに収束するための条件{G(x)}^n=G(a+bx)の解を考えると,n→∞としたときの極値極限分布は,ワイブル分布(指数分布を含む)か2重指数分布(ガンベル分布)の2つのうちどちらかに限られることが証明されています.ワイブル分布は指数分布にしたがう確率変数のベキ乗変換であり,一方,2重指数分布は指数分布にしたがう確率変数の対数変換として導かれる分布です.

 一般の分布について,その極限分布が上にあげた2つの型のいずれかになるというのは非常に有益な結果であり,しかも独立同分布の仮定が多少くずれたとしてもワイブル分布や2重指数分布があてはまることが多いとされています.また,分布の裾が軽い場合の最大値の極限分布としては2重指数分布が現れることも知られていて,たとえば,正規分布f(x)=1/√(2π)exp(-x^2/2)については,y(n)=√(2logn){x(n)-sqr(2logn)}の分布が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

 一様分布f(x)=1/θ 0≦x≦θ F(x)=x/θ

  g(r)=n(n-1)integral(0,θ)(r/θ)^(n-2)(1/θ)^2dx=n(n-1)/θ^nr^(n-2)

 指数分布f(x)=θexp(-θx) F(x)=1-exp(-θx)

  g(r)=n(n-1)integral(0,∞)(exp(-θx)-exp(-θ(x+r)))^(n-2)θexp(-θx)θexp(-θ(x+r))dx

 母集団分布が一様分布や指数分布のときこれは簡単な形になりますが,正規母集団のときには簡単にはなりません.なお,2重指数分布(ガンベル分布)の範囲の分布については第2種変形ベッセル関数を用いて表されることを申し添えておきます.

−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−

(4)標本中央値の分布

 簡単のために,標本の大きさnを奇数とし,n=2m+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/π[tan-1(x-μ)/α]+1/2

より,その標本中央値分布は

  g(x)=(2m+1)!/(m!)2π2^2m{1-(2/πtan-1(x-μ)/α)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})=2√nf(μm)(x-μm)

が漸近的にN(0,1)にしたがうことを示せばよい.

 uの確率密度関数は,

  x=u/2√nf(μm)+μm,dx=du/2√nf(μm)

より,

h(u)du=g(x)dx=g(u/2√nf(μm)+μm)/2√nf(μm)du

h(u)=g(u/2√nf(μm)+μm)/2√nf(μm)

=(2m+1)!/(m!)2*{∫(-∞-u)f(t)dt*∫(u-∞)f(t)dt}^mf(u/2√nf(μm)+μm)/2√nf(μm)

 ここで,n→∞のときの極限を考える.スターリングの法則

  m!=√2πm^(m+1/2)exp(-m)

  (2m+1)!=√2π(2m+1)^(2m+3/2)exp(-2m-1)

を使って簡約化すると

  (2m+1)!/(m!)2→1/√2π*√n*2^n

 また,積分学における第1平均値定理により

  {∫(-∞-u)f(t)dt*∫(u-∞)f(t)dt}^m

=2^(-2m){1-[uf(ξ)]2/n[f(μm)]2}^m→2/2^nexp(-u2/2)

f(u/2√nf(μm)+μm)/2√nf(μm)→1/2√n

したがって,h(u)→1/√(2π)exp(-u^2/2)〜N(0,1)

 さらに,F(q)=pなる統計量qの漸近分布はN(q,p(1-p)/{n[f(μm)]^2})となることを示すこともできます.ここで,p=1/2とおくと中央値の漸近分布,p=1/4とおくと第1四分位数の漸近分布が得られます.また、順序統計量の線形結合(L統計量)の分布に対しても,漸近正規性が成り立つことが証明されています.

 スターリングの公式n!=√2πn^(n+1/2)exp(-n)より

  (2n,n)〜2^(2n)/√(πn)

 母集団分布が正規分布N(μ,σ^2)のとき,

  f(x)=1/√2πσexp{-(x-μ)^2/2σ^2}

ですから,

  μm=μ

  f(μm)=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の正規分布になることも導き出されます.これはα2/(n+2)(π/2)^2と漸近的に等しくなります.したがって,前述の式は簡単な推定方式ながら,かなりよい推定量を与えてくれることがわかります.

===================================