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

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

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

ｆ（ｘ）＝１／π・α／（α^2＋（ｘ−μ）^2）

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

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

∂ｌｎＬ／∂μ＝−２Σ（ｘi−μ）／（（ｘi−μ）^2＋α^2）＝０

∂ｌｎＬ／∂α＝ｎ／α−２αΣ１／（（ｘi−μ）^2＋α^2）＝０

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

ｆ（μ，α）＝∂ｌｎＬ／∂μ＝０

ｇ（μ，α）＝∂ｌｎＬ／∂α＝０

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

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

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

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

ｂ）標本平均値の分散は無限大であるが，標本中央値はμ，標本中央値の分散は有限の値となる

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

　大きさの順に並べ換えた標本（順序統計量）の第１四分位数（２５％点）をＱ1，中央値（５０％点）をＱ2，第３四分位数（７５％点）をＱ3とすると，それぞれμ−α，μ，μ＋αの推定値になりますから，初期値μ0，α0は，分位点を利用してμ0＝Ｑ2または（Ｑ1＋Ｑ3）／２，α0＝（Ｑ3−Ｑ1）／２とします．

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

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

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

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

ｌｎＬ＝-nlogσ-1/σΣ|x−μ|です．

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

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

∂ｌｎＬ／∂σ=-n/σ+1/σ^22Σ|x−μ|=0

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

σの標準誤差は２階偏微分∂2ｌｎＬ／∂σ2より求めます．

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

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

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

ｆ（ｘ）＝１／√２πσｘ・ｅｘｐ（−（ｌｏｇｘ−μ）^2／２σ^2）

　正の値のみをとる確率変数ｘを対数変換してｙ＝ｌｏｇｘとすると，−∞＜ｙ＜∞で，ｙは正規分布Ｎ（μ，σ^2）に従いますから，その密度関数は

ｇ（ｙ）＝１／√２πσ・ｅｘｐ（−（ｙ−μ）^2／２σ^2）で示されます．ここで，∫ｇ（ｙ）ｄｙ＝１，ｄｙ＝ｄｘ／ｘより対数正規分布の密度関数の式が求められます．

【モーメント法】

この分布の母平均，母分散はそれぞれｅｘｐ（μ＋σ^2／２），ｅｘｐ（２（μ＋σ^2／２））（ｅｘｐ（σ^2）−１）で与えられますから，

ｅｘｐ（μ＋σ^2／２）＝ｘ

ｅｘｐ（２（μ＋σ^2／２））（ｅｘｐ（σ^2）−１）＝ｓ^2

を解いて

σ^2＝ｌｎ（（ｘ^2＋ｓ^2）／ｘ^2）

μ＝ｌｎ（ｘ^2／（ｘ^2＋ｓ^2）^1/2）

が求まります．　

【最尤法】

対数の性質から対数正規分布にしたがう変数の積や商は正規分布となります．したがって，ｙi＝ｌｎｘiが平均μ，分散σ^2の正規分布に従うことから，μの最尤推定値はΣｙi／ｎ＝ｌｎn√ｘ1ｘ2・・・ｘn，σ^2の最尤推定値はΣ（ｙi−ｙ）^2／ｎとなります．

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

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

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

ｆ（ｘ）＝１／√２πσ（ｘ−γ）・ｅｘｐ（−（ｌｏｇ（ｘ−γ）−μ）^2／２σ^2）

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

【モーメント法】

ｅｘｐ（μ＋σ^2／２）＋γ＝ｍ1・・・（１）

ｅｘｐ（２（μ＋σ^2／２））（ｅｘｐ（σ^2）−１）＝ｍ2・・・（２）

ｅｘｐ（２（μ＋σ^2／２））（ｅｘｐ（σ^2）−１）^2（ｅｘｐ（σ^2）＋２）＝ｍ3・・・（３）

（３）／（２）よりμを消去してσを求め，（２）に代入してμ，さらに代入してγを求めます．

【最尤法】

　尤度関数は

Ｌ（μ，σ，γ）＝Π１／√２πσ・１／（ｘi−γ）・ｅｘｐ（−（ｘi−γ）^2／２σ^2）

と定義され，∂ｌｎＬ／∂μ＝０，∂ｌｎＬ／∂σ＝０，∂ｌｎＬ／∂γ＝０より

１／σ^2（Σｌｎ（ｘi−γ）／ｎ−μ）＝０，

−ｎ／σ＋１／σ^3（Σｌｎ（ｘi−γ）／ｎ−μ）^2＝０，

Σ１／（ｘi−γ）＋１／σ^2Σ（ｌｎ（ｘi−γ）−μ）／（ｘi−γ）＝０

が得られます．

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

ａ）γの初期値としてデータの最小値より少し小さい値を与える，

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

が考えられますが，著者の経験では，ａ）のほうがｂ）よりもかなり収束安定性が優れていると判断されました．

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

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

ｆ（ｘ）＝１／（αΓ（ｍ））・（ｘ／α）^m-1・ｅｘｐ（−ｘ／α）

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

【モーメント法】

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

平均：ｍα＝ｘ

分散：ｍα^2＝ｓ^2

よりｍ＝ｘ^2／ｓ^2，α＝ｓ^2／ｘ

【最尤法】

　尤度関数は

Ｌ＝Π１／（αΓ（ｍ））・（ｘi／α）^m-1・ｅｘｐ（−ｘi／α）

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

∂ｌｎＬ／∂ｍ＝−ｎ（ｌｎα＋Γ’（ｍ）／Γ（ｍ））＋Σｌｎｘi＝０

∂ｌｎＬ／∂α＝−ｎｍ／α＋Σｘi／α^2＝０

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

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

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

ｆ（ｘ）＝１／（αΓ（ｍ））・（（ｘ−γ）／α）m-1 ・ｅｘｐ（−（ｘ−γ）／α）

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

【モーメント法】

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

平均：γ＋ｍα＝ｍ1

分散：ｍα^2＝ｍ2

平均値まわりの３次モーメント：２ｍα^3＝ｍ3

以上により，

ｍ＝４ｍ2^3／ｍ3^2，α＝ｍ3／２ｍ2，γ＝ｍ1−２ｍ2^2／ｍ3

【最尤法】

　ｍ，α，γは，尤度方程式より，∂ｌ／∂ｍ＝０，∂ｌ／∂α＝０，∂ｌ／∂γ＝０，すなわち

Σｌｎ（ｘi−γ）−ｎｌｎα−ｎΓ’（ｍ）／Γ（ｍ）＝０，

Σ（ｘi−γ）−ｎｍα＝０，

Σ１／（ｘi−γ）−ｎ／｛α（ｍ−１）｝＝０

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

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

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

ｆ（ｘ）＝ｍ／α（ｘ／α）^m-1ｅｘｐ（−（ｘ／α）^m）

　ワイブル分布の尺度母数モデルは上式の関数形で与えられます．この式でｍは形状母数，αは尺度母数と呼ばれます．また，ワイブル分布に従う確率変数を対数変換したものの分布が二重指数分布です．ｙ＝ｌｏｇｘとすると，

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

すなわち，位置母数ｌｏｇα，尺度母数１／ｍの２重指数分布が得られます（ワイブル確率紙はこのことを利用している！）．このことを利用して積率法を適用するほうが簡単なのですが，ここでは馬鹿正直に正攻法でいってみることにします．

【モーメント法】

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

平均：αΓ（１＋１／ｍ）＝ｘ

分散：α^2［Γ（１＋２／ｍ）−｛Γ（１＋１／ｍ）｝^2］＝ｓ^2

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

ＣＶ＝［Γ（１＋２／ｍ）−｛Γ（１＋１／ｍ）｝^2］／Γ（１＋１／ｍ）

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

ｇ（ｍ）＝［Γ（１＋２／ｍ）−｛Γ（１＋１／ｍ）｝^2］／｛Γ（１＋１／ｍ）｝^2−ＣＶ^2＝０

なる式を満たすｍを求めることに帰着します．ＣＶの推定値としては標本変動係数＝ｓ／ｘで代用します．ただし，この標本変動係数を用いる方法はｍの値によっては効率が悪くなります．

【最尤法】

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

ｙ＝Σｘi^mｌｎｘi／Σｘi^m−１／ｍ−Σｌｎｘi／ｎ＝０

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

ｄｙ／ｄｍ＝（Σｘi^mｌｎｘi）^2Σｘi^m−（Σｘi^mｌｎｘi ）^2）／（Σｘi^m）^2＋１／ｍ^2

においてｍの初期値を１としてニュートン・ラプソン法による逐次解法でｍを求め，α^m＝Σｘi^m／ｎよりαを求めます．

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

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

位置=尺度母数モデル

ｆ（ｘ）＝ｍ／α（（ｘ−γ）／α）^m-1ｅｘｐ（−（（ｘ−γ）／α）^m）

【モーメント法】

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

αΓ（１＋１／ｍ）＋γ＝ｍ1・・・（１）

α^2（Γ（１＋２／ｍ）−（Γ（１＋１／ｍ））^2）＝ｍ2・・・（２）

α^3（Γ（１＋３／ｍ）−３Γ（１＋２／ｍ）Γ（１＋１／ｍ）＋２（Γ（１＋１／ｍ））^3）＝ｍ3・・・（３）

ここで（３）／（２）^3/2を計算すればｍのみの関数が得られます．ｍを求めた後，（２）に代入するとαが求まり，（１）に代入するとγが求まります．

【最尤法】

　対数尤度関数は，

ｌｎＬ＝ｎｌｎｍ−ｎｌｎα＋（ｍ−１）Σｌｎ（（ｘ−γ）／α）−Σ（（ｘ−γ）／α）^m

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

ｆ（ｍ，α，γ）＝∂ｌｎＬ／∂ｍ＝０，

ｇ（ｍ，α，γ）＝∂ｌｎＬ／∂α＝０，

ｈ（ｍ，α，γ）＝∂ｌｎＬ／∂γ＝０

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

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

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

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

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

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

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

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

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

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝