■観測データから確率密度関数を求める（その１）
統計的モデリングにおいては，まず，データに何らかの理論分布を想定するという作業が待っています．データの度数分布図を分布関数で規定することができると，データが得られた必然性について確率論的に意味付けすることがが可能になるからです．そして，仮定した数理モデルを検証して理論値が現実の観測値とよく一致しているならば，そのモデルは観察された現象に理論的な説明を与え，データの背後に潜む法則や規則性を発見し，より進んだ結果を予測したりあるいは決定に役立てることができます．
また，２組以上のデータが得られた場合，平均値や分散に差があるか，同分布であるかどうかを問う問題は，最終的にはモデルとして仮定する確率分布の位置母数・尺度母数・形状母数を求める問題に帰着されます．このような現実的要請は多くあり，ここでは未知母数の統計的母数推定法について考えてみることにします．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【１】統計的母数推定法
ｎ個の測定値（ｘ1，ｘ2，・・・，ｘn）に対して，ｍ個の未知のパラメータ（θ1，θ2，・・・，θm）を含む確率密度関数f(x)をあてはめることにします．すなわち，目的は観測データから確率分布の未知パラメータの値を推定することですが，推定には
ａ）パーセント点（分位点）による方法
ｂ）積率推定法（モーメント法）
ｃ）最小χ^2法
ｄ）最尤推定法（最尤法）
などが考案されています．
これらの中で，最小χ^2法と最尤法は区間推定が可能ですが，パーセント点による方法と積率法は点推定のみで区間推定できません．しかも，パーセント点による方法はかなりアバウトな推定法であり，統計的母数推定法の立場からは論外の手法です．
最小χ^2法は一種の重みつき最小２乗法で，期待度数と観測度数から計算されるχ^2値を最小にするようにパラメータの値を定めます．最小χ^2法は一見合理的に思える母数推定法ですが，ヒストグラムは区間の幅の取り方によって見かけが大きく異なってくるため，母数の推定値は区間の作り方に依存して変動するという欠点があります．さらに，最小χ^2法ではデータが階級分割数までに縮約され，階級分割の過程でデータのもっている情報の一部が失われる欠点もあり，この方法をポアソン分布など離散分布なども含めて一般的に適用する際には厄介な問題が生ずることも指摘されています．
そのため，この中でもっとも汎用されている母数推定法は積率法と最尤法ですが，積率法はコンピュータが手軽に利用できなかった頃の計算方法で，現在，そのままの形で利用するには少々古典的であり，最尤法の初期値を求めるのに利用されているにすぎません．
したがって，現在では，積率法でまず点推定し，その値をもとにして最尤法を施してなるべく精密な推定値を誤差も含めて求めるのが，ひとつの標準解法とされています．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【２】最尤法（最大尤度推定法）
最尤法では，標本が与えられたときに事後的にその標本値をもたらした同時確率分布関数（尤度関数）の値を最大にするパラメータの値，すなわち最尤推定値を求めます．同時確率分布関数は互いに独立なｎ個の測定値の確率密度の積
L(θ)=f(x1)f(x2)･･･f(xn)=Πf(xi)
で表され，また，尤度関数の自然対数を対数尤度関数(logarithmic likelihood function)と呼びます．したがって，対数尤度関数は
l(θ)=log(L(θ))=Σlnf(xi)
と定義されます．
対数の単調性から，尤度関数を最大にするパラメータの値と対数尤度関数を最大にするパラメータの値は等しくなります．また，対数尤度関数の値が最大になるパラメータの値は連立方程式
∂l/∂θj=0 （ｊ＝１〜ｍ）
を満たします．この連立方程式を尤度方程式（likelihood equation）といいます．
解析的に母数の最尤推定値を求められる例として，平均μ，分散σ^2の正規分布Ｎ（μ，σ^2）からｎ個の標本（ｘ1，ｘ2，・・・，ｘn）が観測されたという状況を考えて，パラメータμ，σの値を求めてみることにしましょう．
ｘの確率密度関数は
f(x)=1/√2πσexp{-(x-μ)^2/2σ^2)
ですから，尤度関数は
f(x1)f(x2)･･･f(xn)=Πf(xi)
=(2πσ^2)^-n/2exp{-Σ(xi−μ)^2/2σ^2}
ここで両辺の対数をとると
L(μ,σ^2)=Σlnf(xi)
=-2/ln2πσ^2-Σ(xi-μ)^2/2σ^2
になります．
対数尤度を最大にするμとσは
∂L/∂μΣ(xi-μ)/σ^2=0,
∂L/∂σ=-n/σ＋Σ(xi-μ)^2/2σ^3=0
の解としてあたえられますから，したがって，μの最尤推定値はΣxi/n=x（標本平均），σの最尤推定値は√Σ(xi-x)^2/n（標本標準偏差）が得られます．すなわち，母平均の最尤推定値は標本平均，母標準偏差の最尤推定値は標本標準偏差です．
なお，標準偏差σでなく，分散σ^2の最尤推定値を考えると，σ^2を１つの文字とみなして計算すればよいことになり，
∂L/∂σ^2＝−n/2σ^2+Σ(xi-μ)^2/2σ^4=0
これより分散σ^2の最尤推定値はΣ(xi-x)^2/n=s^2（標本分散）となります．つまり，通常使用している標本平均，標本分散，標本標準偏差は最尤推定量となっていることが導かれます．このことから，従来から何気なく使ってきた平均値と分散の理論的背景には，実は最尤推定法があることが理解されるでしょう．
以上のように，尤度さえ表現できてしまえば，あとはその最大化という数値計算の問題に帰着されます．最尤推定値を求めるには対数尤度関数をパラメータに関して微分して０とおきその解を求めますが，正規分布のような例外を除き，推定値を解析的な式の形として得ることは一般には不可能で，専ら数値計算により求めます．１組の非線形方程式の解は，通常，ある初期値を出発点として近似の繰り返し計算によって求めざるをえません．実際には，パラメータの変化分δθjを未知数として，ニュートン・ラプソン法など反復近似法によって解を求めることになります．
１変数非線形方程式の場合のニュートン・ラプソン法の更新公式
Δｘ＝−f(x)/f'(x)
はよく知られていますが，多変数の場合の更新公式は
Σ∂^2L/∂θj∂θkδθk=∂L/∂θj
となります．実際に，L=Σlnf(xi)を代入すると，
∂L/∂θj=Σ1/fi(x)∂fi/∂θj
∂^2L/∂θj∂θk=Σ{1/fi2(x)∂fi/∂θj∂fi/∂θk-1/fi(x) ∂^2fi/∂θj∂θk}
が得られます．
この計算により近似値の修正値が得られますので，この計算を繰り返すことにより真の値に収束させることが期待できます．このように，解析的に解けない問題を近似の繰り返しによって真の値に近づける手法はコンピュータにとって最も得意な分野ですし，アプリケーションプログラムにまかせることにしましょう．
ただし，数値計算上の振る舞いの悪い連続分布関数やパラメータに制約条件の付いた離散分布関数においては単純なニュートン法では解が得られず，特殊な工夫が必要になることを付記しておきます．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【３】最尤推定値の信頼区間
(1)フィッシャーの情報行列
最尤推定値の信頼区間は，フィッシャーの情報行列（Fisher's information matrix）を使って，近似的に求めることができます．フィッシャーの情報行列とは，フィッシャーの情報量の多次元版で，
I(θ)={Ijk}=-{E(∂^2L/∂θj∂θk)}
で定義されるｍｘｍの正方行列です．ここで，Ｉはinformationの略号，E[∂^2L/∂θj∂θk]は対数尤度の２階偏微分：∂^2L/∂θj∂θkの期待値を示します．
本来は，偏微分の積を使って
I(θ)=E[∂logf(x)/∂θj∂logf(x)/∂θk]
で定義されるのですが，
E{∂logf(x)/∂θj]=0,E[∂logf(x)/∂θk]=0
より，２階偏微分行列
I(θ)=-{E(∂^2logf(x)/∂θj∂θk)}
のほうが計算には便利です．
この行列の逆行列I(θ)-1は分散共分散行列であって，その対角要素はパラメータの精度を与える重要な性質をもっていて，パラメータθiの最尤推定量の分散は
Var[θi]=[I(θ)-1]ii=(Δθi)^2
と表されます．分散共分散行列[I(θ)-1]はすべてのパラメータが独立と考えられる場合は対角行列，また，パラメータ間に強い相関がある場合には非対角要素は無視できなくなり対称行列となります．
しかし，一般に期待値（フィッシャーの情報行列）の計算は困難で，期待値の代わりに実際に得られた観測値（ヘシアン）を用いることになります．すなわち，
C(θ)={Cjk}=-{∂^2L/∂θj∂ θk}
をフィッシャーの情報行列の代用とするのです．ｃはcurvature(曲率)の略号で，期待値の代わりに観測値を用いることによって，精度がどのようになるかは，たとえばEfron and Hinckley(1978)を参照されたい．
−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
もう一度，正規分布の場合に戻って，パラメータの標準偏差（標準誤差）Δθを求めてみることにしましょう．フィッシャーの情報行列がどうなるかというと
∂2L/∂μ2=-n/σ^2
∂2L/∂μ∂σ=-2Σ(xi-μ)/σ^3=0
∂2L/∂σ2=-3Σ(xi-μ)^2/σ^4+n/σ^2=-2n/σ^2
ですから，
C={n/σ^2 0 }
[0 2n/σ^2}
したがって，分散共分散行列は
C-1={σ^2/n 0 }
{ 0 σ^2/2n}
対角要素がそのパラメータの標準誤差ですから，
Δμ=σ/√n
Δσ=σ/√2n
となります．
また，正規分布の場合，∂2L/∂μ∂σ=0より，μとσの間には相関は存在しません．
−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
(2)最尤推定量の漸近正規性・漸近有効性
真の母数分布を得るには，それぞれの確率モデルのもとで母数の分布を求める必要がありますが，組合せ数の大きさからいって，このことは非現実的であり，ｎが大きいときの漸近分布を用いることになります．実は，最尤推定量は漸近的に正規分布にしたがう，すなわち，母数θの母分布がたとえ正規分布でなくとも，データ数が多ければ近似的に正規分布からの標本とみなすことができます．これを漸近正規性（asymptotic normality）といいます．
また，不偏推定量のうちで，分散最小の推定量は最小分散不偏推定量と呼ばれ，よい推定量であるとみなされます．したがって，分散の小さい推定量を探してくることが計算統計学の１つの課題となりますが，それでは，分散はどこまで小さくなるのでしょうか？ 実はこれには下限があり，クラーメル・ラオの不等式の下限（Cramer-Rao bound：ＣＲＢ）より母数の誤差はこれ以上小さくはできません．
標本数ｎが大きくなると，最尤推定量の分散は，クラーメル・ラオの下限にするという性質があります．これを漸近有効性（asymptotic efficiency）といいます．すなわち，統計学において母数を推定するのに用いられる最尤法は，漸近的に最小分散不偏推定量と同じ振る舞いをみせる推定量を与えてくれるのです．
−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−−
以上述べたように、最尤推定量は一般に漸近正規性と漸近有効性が成り立ちます．この性質を利用して，最尤推定値の信頼区間（confidence interval）を求めてみましょう．
漸近有効性と漸近正規性を考慮すると，標本数ｎが大きいとき，母数θの100(1-α)%両側信頼区間は，θi±u(α/2)Δθiによって与えられることになります．u(α/2)は正規分布の片側100xα/2%点です．しかし，母数の分布が漸近的に正規分布であるといっても，標本数は常に大きいとは限りません．
母数の分布は厳密に正規分布ではありませんが，幸い，それに近い分布密度になりますから，このような場合，母数の分布としてｔ分布を想定するのが常套手段です．ｔ分布は自由度ν→∞の場合は正規分布を表わしますから，正規分布を一般化したものと考えることができます．
したがって，標本数が小さい場合も考慮して，θiの100(1-α)%両側信頼区間をθi±t(α/2,n-m)Δθiとします．ここで，t(α/2,n-m)は自由度n-mのｔ分布の片側100xα/2%点で，u(α/2)より少し広めに信頼区間を設定したことになります．
自由度をn-mとしたのは漸近不偏性への対処と考えられますが、厳密正規性でなく，漸近正規性しか仮定できない場合は，この方法が一般的に行われています．
しかし、さらに厳密に考えるとが、漸近正規性を有する分布はｔ分布ばかりではありません。一般的ではありませんが、最尤推定量であることを考えると、本来はｔ’分布を用いるべきと思われます。ｔ’分布はｔ分布よりさらに幅の広い分布ですから，θi±t'(α/2,n-m)Δθiとおくと、さらに広めに信頼区間を設定したことになります．
t'(α/2,n-m)値はt(α/2,n-m)値から簡単に換算できます．少し控えめすぎるかもしれませんが，期待値の代わりに観測値を用いているわけですから，ここではこの結果を採用することにします．
繰り返しますが，最尤法で求めた母数θjの分布形は，厳密にいうと正規分布になるとは限りません．ですから，各母集団分布形ごとに母数分布を用意する必要が生じますが，この組合せの数は極めて多数となり，その分布も特殊なものになります．
母数の正確な分布を求めることはｎを適当な大きさにしてもたいそう労力を要することであり，精密な式は少数の例外を除いて求められそうにありません．したがって，信頼区間についてはっきりした議論をしようと思ったら，合わせる関数形ごとにシミュレーションで母数の分布を調べたほうが安全ということになります．
しかし，厳密に正規分布ではないもののそれに近い分布に対しては，それが正規分布するとみなして近似的に求めることができるというのが，ここでの議論です．正規近似の結果，精密度の概念はいくらかぼやけてきますが，実際問題として，母数の分布が極度に非正規であることは少なく，少し広めに信頼区間を設定すればよいわけで，一部は経験的・習慣的考えに基づき，一部は理論的な考えに基づいて議論を展開しました．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
【４】積率法（モーメント法）
一般的にいって，尤度方程式は非線形連立方程式になります．非線形方程式を解く場合，計算がうまく収束するかどうかが初期値の選び方によって決まることが多いため，初期値の見積もりは重要な課題です．幸い，分布関数の母数推定では，積率法を用いて初期値を比較的簡単に求めることができます．
積率法（モーメント法）では，確率分布の理論モーメント（μ1 ，μ2 ，μ3 ，・・・）と標本モーメント（ｍ1 ，ｍ2 ，ｍ3 ，・・・）がそれぞれ等しいとおいた方程式（μk=mk）を確率分布のパラメータの数だけ作ればパラメータの値が推定でき関数形は確定します．つまり，２母数確率分布では２つの方程式，３母数確率分布では３つの方程式が必要になります．ここで，実測値であるｋ次標本モーメントｍkは
ｍk=Σ(x-x)k/nによって計算することができます．
ｍ1=Σx/n（標本平均x)
ｍ2=Σ(x-x)^2/n （標本分散s^2）
ｍ3=Σ(x-x)^3/n
ｍ4=Σ(x-x)^4/n
正規分布の場合について説明すると，分布関数の未知パラメータはμ1＝μ，μ2＝σ^2の２母数ですから，おのおのに対してその推定量であるｍ1（標本平均ｘ）とｍ2（標本分散ｓ^2）を代入すれば未知母数が決定できます．
μ＝ｘ（標本平均）
σ^2＝ｓ^2（標本分散）
正規分布では積率法と最尤法は同じ推定量を与えることがわかりましたが，積率法ではより容易に推定量を求められるため，自然の要求にあった推定法といえます．しかし，積率法では点推定しかできず，精度保証つきという面からいって区間推定できる最尤法のほうが格段に優れた強みをもっています．
＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝
