■観測データから確率密度関数を求める(その2)

【1】非正規分布における積率法と最尤法

 積率法,最尤法は確率密度関数ごとに個別に対応していかなければなりません.代表的な分布をいくつか選んで,その方法を記述してみることにしましょう.

(1)コーシー分布の母数推定

f(x)=1/π・α/(α^2+(x−μ)^2)

 コーシー分布はt分布において自由度1としたものであり,平均値は定まらず分散が無限大になる厄介な分布です(平均や分散をもたない確率分布!).したがって,コーシー分布ではモーメント法を適用して分布のパラメータを求めることはできず,パラメータの推定には専ら最尤法を適用することになります.

 対数尤度関数の偏微分を0とおくと

∂lnL/∂μ=−2Σ(xi−μ)/((xi−μ)^2+α^2)=0

∂lnL/∂α=n/α−2αΣ1/((xi−μ)^2+α^2)=0

になりますから,最尤法では

f(μ,α)=∂lnL/∂μ=0

g(μ,α)=∂lnL/∂α=0

とおいて漸近連立方程式の逐次解法によりμとαを計算します.

 また,コーシー分布に従う変数については,中心極限定理「母集団分布が正規分布でなくても標本が大きくなると標本平均値の分布は次第に正規分布に近づく」が成り立ちません.一方,標本中央値の分散は標本の大きさnを大きくすると小さくなることが示されます.

 母数μ,αの推定に際して,初期値μ0,α0のとり方をやみくもに行なうと収束が非常に悪くなりますが,

a)コーシー分布では平均値は定まらなくても中央値,最頻値はμである

b)標本平均値の分散は無限大であるが,標本中央値はμ,標本中央値の分散は有限の値となる

ことから,コーシー分布に対して最尤法を適用する場合は,標本平均値でなく標本中央値を初期値として与えパラメータ推定を行ないます.

 大きさの順に並べ換えた標本(順序統計量)の第1四分位数(25%点)をQ1,中央値(50%点)をQ2,第3四分位数(75%点)をQ3とすると,それぞれμ−α,μ,μ+αの推定値になりますから,初期値μ0,α0は,分位点を利用してμ0=Q2または(Q1+Q3)/2,α0=(Q3−Q1)/2とします.

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

(2)両側指数分布(ラプラス分布)

f(x)=1/2σexp(-|x−μ|/σ)

 両側指数分布の対数尤度関数は,定数部分を省略すると

lnL=-nlogσ-1/σΣ|x−μ|です.

 この関数の偏微分は絶対値を含むことになり,正攻法で解こうとするとなかなか厄介です.そこで,対数尤度関数の偏微分を0とおくのではなく,直接,対数尤度関数を最大にすることを考えます.もしσが既知ならば,対数尤度関数を最大にするにはΣ|x−μ|を最小にしなければなりませんから,μの最尤推定量は標本中央値になることがわかります.分散を最小にする統計量が平均であるのに対し,絶対偏差を最小にする統計量はメジアンであるからです.

 ラプラス分布の位置母数の最尤推定量は標本中央値であり,μが求まったところでσで偏微分すると

∂lnL/∂σ=-n/σ+1/σ^22Σ|x−μ|=0

したがって,σ=1/nΣ|x−μ|

σの標準誤差は2階偏微分∂2lnL/∂σ2より求めます.

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

(3)2母数対数正規分布の母数推定

 対数変換したものが正規分布となる分布が対数正規分布で,その確率分布は次式で表現されます.

f(x)=1/√2πσx・exp(−(logx−μ)^2/2σ^2)

 正の値のみをとる確率変数xを対数変換してy=logxとすると,−∞<y<∞で,yは正規分布N(μ,σ^2)に従いますから,その密度関数は

g(y)=1/√2πσ・exp(−(y−μ)^2/2σ^2)で示されます.ここで,∫g(y)dy=1,dy=dx/xより対数正規分布の密度関数の式が求められます.

【モーメント法】

この分布の母平均,母分散はそれぞれexp(μ+σ^2/2),exp(2(μ+σ^2/2))(exp(σ^2)−1)で与えられますから,

exp(μ+σ^2/2)=x

exp(2(μ+σ^2/2))(exp(σ^2)−1)=s^2

を解いて

σ^2=ln((x^2+s^2)/x^2)

μ=ln(x^2/(x^2+s^2)^1/2)

が求まります. 

【最尤法】

対数の性質から対数正規分布にしたがう変数の積や商は正規分布となります.したがって,yi=lnxiが平均μ,分散σ^2の正規分布に従うことから,μの最尤推定値はΣyi/n=lnn√x1x2・・・xn,σ^2の最尤推定値はΣ(yi−y)^2/nとなります.

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

(4)3母数対数正規分布の母数推定

 対数正規分布においてデータxに下限がある場合には,位置母数をγとして,x−γをxの代わりに用いて理論分布をあてはめることが必要になってきます.そこで,位置母数γを追加した3母数対数正規分布

f(x)=1/√2πσ(x−γ)・exp(−(log(x−γ)−μ)^2/2σ^2)

の未知母数μ,σ,γを推定する手順を記しておきます.

【モーメント法】

exp(μ+σ^2/2)+γ=m1・・・(1)

exp(2(μ+σ^2/2))(exp(σ^2)−1)=m2・・・(2)

exp(2(μ+σ^2/2))(exp(σ^2)−1)^2(exp(σ^2)+2)=m3・・・(3)

(3)/(2)よりμを消去してσを求め,(2)に代入してμ,さらに代入してγを求めます.

【最尤法】

 尤度関数は

L(μ,σ,γ)=Π1/√2πσ・1/(xi−γ)・exp(−(xi−γ)^2/2σ^2)

と定義され,∂lnL/∂μ=0,∂lnL/∂σ=0,∂lnL/∂γ=0より

1/σ^2(Σln(xi−γ)/n−μ)=0,

−n/σ+1/σ^3(Σln(xi−γ)/n−μ)^2=0,

Σ1/(xi−γ)+1/σ^2Σ(ln(xi−γ)−μ)/(xi−γ)=0

が得られます.

 ニュートン・ラプソン法により,この方程式を解いて,μ,σ,γの最尤推定量を求めます.その際,γの初期値の求め方には2つの方法

a)γの初期値としてデータの最小値より少し小さい値を与える,

b)モーメント法で求めたγを初期値とする

が考えられますが,著者の経験では,a)のほうがb)よりもかなり収束安定性が優れていると判断されました.

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

(5)2母数ガンマ分布の母数推定

f(x)=1/(αΓ(m))・(x/α)^m-1・exp(−x/α)

 ガンマ分布は2個のパラメータm,αを用いて,上式で定義される分布です(尺度母数モデル).αはスケールの取り方に関係するパラメータ(尺度母数)で,厳密な意味では曲線の形に影響を及ぼすものではなく,単に縦横の座標を伸縮するのと同じ意味をもつにすぎません.mのほうは曲線の本質的な形を決定するもの(形状母数)で,αに比して重要な意義をもっています.

【モーメント法】

ガンマ分布では,モ−メント法で容易に推定量を求めることができます.

平均:mα=x

分散:mα^2=s^2

よりm=x^2/s^2,α=s^2/x

【最尤法】

 尤度関数は

L=Π1/(αΓ(m))・(xi/α)^m-1・exp(−xi/α)

となるので,尤度方程式は

∂lnL/∂m=−n(lnα+Γ’(m)/Γ(m))+Σlnxi=0

∂lnL/∂α=−nm/α+Σxi/α^2=0

で与えられます.モーメント法で求めたm,αを初期値として母数を推定するわけですが,ニュートン・ラプソン法など非線形連立方程式の逐次解法を用いることになります.

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

(6)3母数ガンマ分布の母数推定

f(x)=1/(αΓ(m))・((x−γ)/α)m-1 ・exp(−(x−γ)/α)

 2母数ガンマ分布は尺度母数モデルだったわけですが,これに位置母数をつけて位置=尺度母数モデルにしたものがこの分布です.

【モーメント法】

 3母数ガンマ分布とは,位置母数つきのガンマ分布を指しますが,前節のように,モ−メント法を適用させます.

平均:γ+mα=m1

分散:mα^2=m2

平均値まわりの3次モーメント:2mα^3=m3

以上により,

m=4m2^3/m3^2,α=m3/2m2,γ=m1−2m2^2/m3

【最尤法】

 m,α,γは,尤度方程式より,∂l/∂m=0,∂l/∂α=0,∂l/∂γ=0,すなわち

Σln(xi−γ)−nlnα−nΓ’(m)/Γ(m)=0,

Σ(xi−γ)−nmα=0,

Σ1/(xi−γ)−n/{α(m−1)}=0

の解として求められます.しかしながら,mが1に近い値の場合にはこの非線形連立方程式の解は不安定になります.

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

(8)2母数ワイブル分布の母数推定

f(x)=m/α(x/α)^m-1exp(−(x/α)^m)

 ワイブル分布の尺度母数モデルは上式の関数形で与えられます.この式でmは形状母数,αは尺度母数と呼ばれます.また,ワイブル分布に従う確率変数を対数変換したものの分布が二重指数分布です.y=logxとすると,

g(y)=mexp[m(y-logα)-exp{m(y-logα)}]

すなわち,位置母数logα,尺度母数1/mの2重指数分布が得られます(ワイブル確率紙はこのことを利用している!).このことを利用して積率法を適用するほうが簡単なのですが,ここでは馬鹿正直に正攻法でいってみることにします.

【モーメント法】

 モーメント法による2母数ワイブル分布のパラメータ推定では

平均:αΓ(1+1/m)=x

分散:α^2[Γ(1+2/m)−{Γ(1+1/m)}^2]=s^2

なる形状母数mと尺度母数αを求めればよいのですが,母変動係数

CV=[Γ(1+2/m)−{Γ(1+1/m)}^2]/Γ(1+1/m)

がmのみに依存することを用いると,mの推定は,

g(m)=[Γ(1+2/m)−{Γ(1+1/m)}^2]/{Γ(1+1/m)}^2−CV^2=0

なる式を満たすmを求めることに帰着します.CVの推定値としては標本変動係数=s/xで代用します.ただし,この標本変動係数を用いる方法はmの値によっては効率が悪くなります.

【最尤法】

 次に,ワイブル分布のパラメータ推定に最尤法を適用してみます.尤度関数を偏微分して0と置いた式を整理するとmは

y=Σxi^mlnxi/Σxi^m−1/m−Σlnxi/n=0

の解として求められます.

dy/dm=(Σxi^mlnxi)^2Σxi^m−(Σxi^mlnxi )^2)/(Σxi^m)^2+1/m^2

においてmの初期値を1としてニュートン・ラプソン法による逐次解法でmを求め,α^m=Σxi^m/nよりαを求めます.

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

(9)3母数ワイブル分布の母数推定

位置=尺度母数モデル

f(x)=m/α((x−γ)/α)^m-1exp(−((x−γ)/α)^m)

【モーメント法】

 3母数ワイブル分布では,平均値と分散だけでなく与えられたデータの歪度を用いてパラメータの値を求めることにします.

αΓ(1+1/m)+γ=m1・・・(1)

α^2(Γ(1+2/m)−(Γ(1+1/m))^2)=m2・・・(2)

α^3(Γ(1+3/m)−3Γ(1+2/m)Γ(1+1/m)+2(Γ(1+1/m))^3)=m3・・・(3)

ここで(3)/(2)^3/2を計算すればmのみの関数が得られます.mを求めた後,(2)に代入するとαが求まり,(1)に代入するとγが求まります.

【最尤法】

 対数尤度関数は,

lnL=nlnm−nlnα+(m−1)Σln((x−γ)/α)−Σ((x−γ)/α)^m

で,尤度関数を最大にする必要条件は

f(m,α,γ)=∂lnL/∂m=0,

g(m,α,γ)=∂lnL/∂α=0,

h(m,α,γ)=∂lnL/∂γ=0

なる非線形連立方程式で与えられますから,漸近連立方程式の逐次解法であるニュートン・ラプソン法を用いてm,α,γを求めます.

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

【2】モーメント法と最尤法の比較

 統計学は,観察データに基づく記述統計学と,実験に基づく小標本により母集団のパラメータを推測する推測統計学に二分されます.前者は進化論で有名なダーウィン(自然淘汰説)に礎をおき,カール・ピアソンにより大成されたものです.また,後者は遺伝の法則を導いたメンデル(遺伝学説)の流れをくみ,フィッシャーにより大成されました.

 統計学を応用数学の一部門の水準まで高めること,統計学の応用領域を拡大創建することが,二人の終生の念願でありました.ところが,ピアソンとフィッシャーは両者とも強い個性の持ち主であったようで,生涯に何度も有名な論争をしています.たとえば,観察データに理論分布を適合させる手段として積率法を導入したピアソンに対して,フィッシャーは,標本から母集団のパラメータを推定する際に,モーメント法(積率法)による推定量はよい特性をもっていないと批判し,モーメント法に代わるものとして最尤推定法(最尤法)を提唱しています.

 現在では,モーメント法で求めたパラメータを初期値として最尤法を用いてなるべく精密な推定値を誤差も含めて求めておくのが,ひとつの標準解法になっています.

 積率法と最尤法の性能比較を行ってみましょう.フィッシャーの指摘のごとく,ピアソンのモーメント法では高々数次までのモーメントの一致性が保証されるのみで高次のモーメントでは不偏性はおろか一致性すら成立しないため,その推定量は満足すべき特性をもっていません.また,これまでの議論により,最尤法では,分散共分散行列の対角要素からパラメータの誤差を同時推定することができます.したがって,最尤推定法の一つの型の中で議論できるという強みをもっている最尤法に比較すると,モーメント法では推定効率が低くなり,母数の点推定のみで区間推定することはできません.さらに,積率法はコーシー分布のように積率が存在しない場合には使えないという欠点もあります.

 しかし,最尤法において尤度方程式を解く方法は非常に煩雑であり計算機の助けなしにはほとんど不可能です.また,分布によっては正則条件が満たされず適用できませんし,解も不安定になり,大きな計算誤差を伴うかあるいは計算不能になってしまうため,簡単には未知パラメータを求められなくなることがあります.そのような場合のパラメータ推定には,緊急避難的にモーメント法が使われます.

 最尤法は積率法ほどには直感に訴える意味を持ち合わせておらず,喩えていうならば,前者は論理を優先し面白味に欠けるが着実である方法なのに対し,後者のほうが直感的で人の心にストレートに訴えるものがあるが往々にして失敗する方法というわけです.このように,確率分布のパラメータ推定に最良かつ万能の方法は存在しません.そのため,現実的にはモーメント法と最尤法が相互補完的に用いられているのです.

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