■分布の混合について

 確率変数xが標準正規分布

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

にしたがうとき,y=x^2の分布はどうなるのでしょうか?

  dy=2xdx

また,y=x^2の逆関数は単調増加関数ではなく,2価関数であることに注意すると

  g(y)dy=2f(x)dx

  g(y)2xdx=2/√2πexp{-x^2/2}dx

より,

  g(y)=1/√2πy^(-1/2)exp{-y/2}

が得られます.この分布は自由度1のχ^2分布と呼ばれるものです.

 このように確率変数を変数変換した関数の分布は簡単に求めることができます.今回のコラムでは確率変数の変数変換ではなく,似て非なる問題「分布の混合」について取り上げてみたいと思います.

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

【1】分布の混合

 分布の混合(compoundまたはmixture of probability distribution)という用語は2つの意味で使われています.

 1つは,複数の独立な確率分布

  f1(x),f2(x),・・・,fk(x)

がある比率

  w1,w2,・・・,wk  (Σwi=1,wi>0)

をもって混ぜ合わさっている場合であり,

  f(x)=w1f1(x)+w2f2(x)+・・・+wkfk(x)=Σwifi(x)

が求める確率密度関数となります.これはいわば単純な加重和であって,混合される各分布は成分ということになります.

 たとえば,混合正規分布

  F(x)=(1-ε)Φ(x-μ/σ)+εΦ(x-μ/3σ)

  f(x)=(1-ε)/σφ(x-μ/σ)+ε/3σφ(x-μ/3σ)

は2つの成分よりなり,(1-ε)の確率でN(μ,σ^2)に従い,εの確率でN(μ,9σ^2)に従います.

 また,合計1になる正の重み自体が1つの確率分布(g)に従うと考えることができます.例をあげると,非心χ^2分布はχ^2分布(f)とポアソン分布(g)の混合となっています.

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

 もう1つは,確率密度関数のパラメータが確率変数となっている確率分布関数を指す場合です.母数θを含む確率密度関数f(x)をパラメータまで含めてf(x,θ)と表します.f(x,θ)の母数θがまた確率密度関数g(θ)に従うとすると

  h(x)=∫f(x,θ)g(θ)dθ

によって新しい確率密度関数h(x)が定義されます.

 また,混合分布の累積分布関数は

  H(x)=∫h(x)dx=∫F(x,θ)g(θ)dθ=∫F(x,θ)dG(θ)

となります.

 これでθに関する重み付き平均を考えたことになりますが,位置母数モデルf(x,θ)=f(x-θ)のとき,hはfとgとの畳み込みそのものです.たとえば,X線,γ線などの電磁波はそれぞれの線スペクトルに固有の幅と分布をもっていて,光の線スペクトルのようなコーシー分布(f):

  f(x)=1/π・β/(β^2+(x-θ)^2)

を示すものを分光器で測定したとすると,分光器には固有の分解能があり,それは正規分布(g)

  g(θ)=1/√2πσexp(-θ^2/2σ^2)

で近似できることが多いわけですから,測定したスペクトルの分布hはコーシー分布と正規分布を合成したものになります.すなわち,フォークト関数(h)はコーシー分布と正規分布の位置混合になっているというわけです.

 一方,尺度母数モデルf(x,θ)=f(x/θ)/θのとき,hはfとgの尺度混合であるといいます.g(θ)としてはガンマ分布(カイ2乗分布)が用いられることが多く,たとえば,ポアソン分布(f)をガンマ分布(g)で混合すれば負の2項分布(h),指数分布を(f)をガンマ分布(g)で混合すればパレート分布(h)が得られます.t分布も正規分布(f)のσ^2がガンマ分布(g)にしたがうと仮定,すなわちガンマ分布で尺度混合した混合分布です.

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

[補]畳み込み(convolution)

 ある生体システム(g)にあるパルス(f)を入力した場合,出力としてh(t)が観察された場合,たとえば,薬物経口投与後の血中濃度h(t)は,静注後の血中濃度fと吸収速度gの関数として表現することができます.このようにfとgを合成した関数は畳み込みと呼ばれ,一般に,分布のたたみ込みは

  p*q(x)=∫(-∞,∞)p(x-t)q(t)dt

として定義されます.この式は,薬物体内動態評価のみならず,量子力学,数理生物学,赤外線分析,X線回折,クロマトや反応工学および多くの物理現象にわたって現れる普遍的関係式です.

 ガウス関数(Gaussian,正規分布):

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

とローレンツ関数(Lorentzian,コーシー分布):

  g(x)=1/π・β/(β^2+(x-α)^2)

の畳み込み(convolution)は

  h(z)=∫(-∞,∞)f(z-y)g(y)dy

    =∫(-∞,∞)g(z-x)f(x)dx

    =1/πσ∫(0,∞)exp(-t^2/2-βt/σ)cos{(α+μ-z)/σt}dt

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

【2】正規母集団からの標本分散の分布

 これ以降は後者の場合について考えてみますが,「母集団分布が正規分布であるとき,標本分散と母分散の比はχ^2分布にしたがう」すなわち

  (n-1)u^2/σ^2〜χ^2(n-1)

  ns^2/σ^2〜χ^2(n-1)

という統計的性質を使って,標本分散(の平方根)の分布がどうなるかを考えてみることにします.

 まず最初にuの分布を求めてみることにしましょう.自由度n-1のχ^2分布は

  f(x)dx=1/2^((n-1)/2)Γ((n-1)/2)x^((n-3)/2)exp(-x/2)dx

で表されますから,x=(n-1)u^2/σ^2で変数変換すると,uの分布:

  g(u)du=1/2^((n-1)/2)Γ((n-1)/2)((n-1)u^2/σ^2)^((n-3)/2)exp(-(n-1)u^2/2σ^2)2(n-1)u/σ^2du

が得られます.

 次に,標本分散の期待値と分散を求めてみましょう.

  E[u^k]=∫u^kg(u)du

において,(n-1)u^2/2σ^2=tとおいて整理すると,

  E[u^k]=σ^kΓ((n-1)/2+k/2)/Γ((n-1)/2)(2/(n-1))^(k/2)

が得られます.

 この式はシュワール(Shewhart)の公式と呼ばれます.ガンマ関数の計算になれていないと,式の誘導は難しいのでおしつけがましく結果だけを書いておきますが,ここではとりあえず,信じるものは救われる.ホレ信じなさい.というわけですが,天下り式で我慢できないかたはガンマ関数の記述を参照しながら誘導を試みられたい.

 ともあれ,これより,k=1の場合,

  E[u]=σΓ(n/2)/Γ((n-1)/2)(2/(n-1))^(1/2)

したがって,

  E[uΓ((n-1)/2)/Γ(n/2)((n-1)/2)^(1/2)]=σ

となり,uΓ((n-1)/2)/Γ(n/2)((n-1)/2)^(1/2)がσの不偏推定値です.

 同様にして

  k=2ではE[u^2]=σ^2(σ^2の不偏推定値はu^2)

  k=3ではE[u^3Γ((n-1)/2)/Γ(n/2+1)((n-1)/2)^(1/2)]=σ^3

  k=4ではE[(n-1)/(n+1)u^4]=σ^4

が成立します.不偏推定値において,kが偶数のときは簡単になって

  E[u^k(n-1)^(k/2)/(n+k-3)(n+k-5)・・・(n-1)]=σ^k

が得られます.

 なお,正規分布の場合,母標準偏差の不偏推定量は

  uΓ((n-1)/2)/Γ(n/2)((n-1)/2)^(1/2)

のようになりますが,一般の分布の場合,母標準偏差の不偏推定量は見いだされていません.

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

 標本不偏分散u^2の期待値が母分散σ^2に一致することが確かめられましたが,上のことを利用すると簡単に標本不偏分散の分散が求まります.

  E[u^2]=σ^2

  E[u^4]=σ^4(n+1)/(n-1)

  V[u^2]=E[(u^2-E[u^2])^2]

     =E[u^4]-E[u^2]^2=σ^4(n+1)/(n-1)-σ^4=2σ^4/(n-1)

 一方,標本分散s^2の平方根sは最尤標準偏差ですが,不偏分散u^2の平方根uは不偏標準偏差にはなりません.母標準偏差σの不偏推定量,すなわち,E[d]=σなるdは,

  d=√n/2Γ{(n-1)/2}/Γ(n/2)s=√(n-1)/2Γ{(n-1)/2}/Γ(n/2)u

で与えられます.ここで,sやuの係数は1より大きい値を示すところから,  d>u>s

になります.

 ついでに,σ^4の不偏推定値としてはu^4を使えばよいのだろうかという問題の解答も掲げておきます.結論だけかいておくと

  n^2/(n^2-1)s^4=(n-1)/(n+1)u^4=1/(n^2-1){Σ(x-x)^2}^2

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

【3】スチューデントのt分布

 測定では一般に母分散σ^2を知ることはできません.そのかわり,何回かの測定のばらつきからnで割った標本分散s^2とn−1で割った標本不偏分散u^2が求められます.

  s^2=Σ(x-x)^2/n

  u^2=Σ(x-x)^2/(n-1)

 スチューデントやウェルチのt検定では母分散σ^2の代わりにその推定値を使うことになりますが,この推定量としては標本分散s^2と標本不偏分散u^2が考えられます.nが大きいとき標本分散s^2と標本不偏分散u^2の差は小さくなりますから,あまり問題にはなりませんが,nが小さいときにはその違いは案外大きくなります.

 前節で標本分散s^2は母分散σ^2の最尤推定量,標本不偏分散u^2は不偏推定量になっていることを喚起しておきましたが,それでは母分散σ^2の推定量としてどちらが好ましいのでしょうか?

 標本分散s^2はσ^2の最尤推定量であるだけでなく,その平方根もまたσの最尤推定量です.一方,標本不偏分散u^2の平方根uはσの不偏推定量にはなっていません.標本分散の平方根を用いると最尤推定量の一つの枠の中で議論できるという強みをもっているため,その意味で,u^2でなくs^2を使うべきだと主張する人もいます.

 結論からいうと,母分散σ^2の代わりにはその不偏分散u^2が使われます.(1)「nで割る標本分散」は本来の分散σ^2より小さい方に偏った答えを与える傾向があり,とくに,実験データのようにデータが少数の場合,このような偏りを補正するために「n−1で割る不偏分散」を使うべき

(2)母標準偏差の不偏推定量,

  √n/2Γ{(n-1)/2}/Γ(n/2)s=√(n-1)/2Γ{(n-1)/2}/Γ(n/2)u

の定数はt分布の確率密度関数の定数であるから

というのがもっともらしい理由としてあげられますが,これらは二義的な理由にすぎません.

 t分布とのつながりを考えるとき,不偏分散(の平方根)でないといけないというのが第一義的な理由です.この節ではこのことを証明したいと思います.

(証明)

 n個の観測値の標本平均と母平均の差(距離)を不偏標本分散の平方根で割った統計量t=/√n(x−μ)/uの分布が自由度n−1のt分布に従うことは,ゴセットが最初に発見しフィッシャーが厳密に証明したことは歴史的事実として有名です.フィッシャーは統計量tの分布をn次元ユークリッド空間を使って導きましたが,オリジナルの導出法は多重積分が必要となりわずらわしいので,再び,「母集団分布が正規分布であるとき,標本不偏分散と母分散の比はχ^2分布にしたがう」という統計的性質を使って求めてみることにします.

 前節と同様にして,σの分布は

  g(σ)dσ=-1/2^((n-1)/2)Γ((n-1)/2)((n-1)u^2/σ^2)^((n-3)/2)exp(-(n-1)u^2/2σ^2)2(n-1)u^2/σ^3dσ

で与えられます.

 標本平均xの分布はN(μ,σ^2/n)すなわち

  f(x)=√(n/2π)/σexp(-n(x-μ)^2/2σ^2)

であり,σの分布も上式で与えられますから,分布の混合によって

  h(x)=∫(0,∞)f(x,σ)g(σ)dσ

 =√(n/2π(n-1))Γ(n/2)/Γ((n-1)/2)u^(n-1){n/(n-1)(x-μ)2+u2}^(-n/2)

 ここで,√n(x-μ)/u=tと変数変換すると

  h(μ)dμ=1/√π(n-1)Γ(n/2)/Γ((n-1)/2){1+t^2/(n-1)}^(-n/2)dt

より

  h(t)=1/√π(n-1)Γ(n/2)/Γ((n-1)/2){1+t^2/(n-1)}^(-n/2)

 このようにして母集団の分散によらない分布が構成されます.この分布は自由度n−1のt分布にほかなりません.これより,「xがN(0,1),yが自由度nのχ^2分布に従うとき,t=x/√(y/n)は自由度nのt分布に従う.」という統計的性質が得られます.また,これによりt分布は正規分布よりも裾の重い混合正規分布の1つと考えることができます.

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

【4】ワイブル分布とパレート分布の混合

 前節ではt分布が正規分布とガンマ分布の尺度混合であることを示しましたが,ここではワイブル分布とパレート分布の尺度混合について考えることにします.

 コラム「ポアソン配置とワイブル分布(雑然か整然か)」では,n次元空間中に無作為に配置(ポアソン配置)した点の集団があるとき,ある点からその最も近い点(最近接点:nearest neighbor)に至るまでの距離の分布を考えました.

 任意のn次元における最近接点間の距離の分布は,ポアソン配置した点の集団の密度をδとして,ワイブル分布

  f(x)=sn-1δx^(n-1)exp(-vnδx^n)

      =nvnδx^(n-1)exp(-vnδx^n)

  F(x)=1−exp(-vnδx^n)

になることが誘導されます.ただし,

  vn=π^(n/2)/Γ(n/2+1),v1=2,v2=π,v3=4π/3

 ここで,局所的にポアソン配置と仮定できる点の集団の密度δがパレート分布

  δ(r)=k/r^a

に従って連続的に変化する場合を考えます.すなわち,局所的な分布はポアソン配置とみなせるが大域的な密度は双曲線様のパワー曲線になるというものです.

 分布の混合により

  h(x)=∫(r1,r2)f(x,δ(r))δ(r)dr

 =nk/{a(vnx^nk)^(1-1/a)}{δ(2-1/a,vnx^nk/r1^a)-δ(2-1/a,vnx^nk/r2^a)}

 ここでδ(c,x)は不完全ガンマ関数を表しますが,合流型超幾何関数1F1を用いて

  δ(c,x)=x^c/c1F1(1,1+c,-x)=1-x/(1+c)+x^2/(1+c)(2+c)-x^3/(1+c)(2+c)(3+c)+・・・

あるいは

  δ(c,x)=x^c/c・exp(-x)1F1(1,1+c,x)

で具体的な値を求めることができます.

 r1→0,r2→∞では

  h(x)=nk/{a(vnx^nk)^(1-1/a)}Γ(2-1/a)

となって,最近接距離はパレート分布に従うことがわかります.また,r1→0,r2→rとすると

  h(x)=nk/{a(vnx^nk)^(1-1/a)}{Γ(2-1/a)-δ(2-1/a,vnx^nk/r^a)}

=nk/{ay^(1-1/a)r^(a-1)}{Γ(2-1/a)-1+y/(3-1/a)-y^2/(4-1/a)/(5-1/a)+y^3/(3-1/a)(4-1/a)(5-1/a)-・・・} y=vnx^nk/r^a

 なお,

  h(x)=∫(r1,r2)f(x,δ(r))δ(r)dr

において,δ(r)=δ(定数)のとき

  h(x)=∫(r1,r2)f(x,δ)δdr=(r2-r1)f(x,δ)δ

となって,混合分布はワイブル分布になることを注意しておきます.

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