最尤法とベイズ的方法
２項分布の区間推定
統計学的分析法における代表的なものが最尤法（Maximum Likelihood Estimation, MLE）であり、最近ではベイズ的方法（Bayesian Data Analysis）が注目されている。２項分布を例として、これら２つのプログラム例をVisual C++で用意した。まず、最尤法による分析のプログラムを説明する。最尤法については、「岡本安晴『統計学を学ぶための数学入門[下]』2009、培風館」の「9.12節 確率モデルによるデータ分析」における「最尤法」、「岡本安晴『データ分析のための統計学入門』2009、おうふう」の「6.3節 最尤法」、「岡本安晴『計量心理学』2006、培風館」の「付録Ｆ 尤度関数・最尤法・尤度比・AIC・フィッシャー情報行列」で説明している。最尤法では母数（パラメタ）の推定値を尤度関数を最大にするものとして与える。２項分布の場合の尤度関数は
で与えられるので、母数の推定値は
となる（岡本安晴「統計学を学ぶための数学入門[下]」、Pp. 152-153参照）。
95％信頼区間はフィッシャー情報行列を用いると
となる（岡本安晴「データ分析のための統計学入門」、Pp. 155-160参照）。ここで、
である。
以上のことから、最尤法によって母数の推定値のその95％信頼区間を求めるプログラムmleは以下のようになる。
int N = int::Parse(textBoxN->Text);
int k = int::Parse(textBoxK->Text);
richTextBox1->AppendText("\nN = " + N.ToString() + "
k = " + k.ToString());
double
p_hat = double(k) / N;
double
sigma2 = p_hat * (1.0 - p_hat);
double
hfwidth95 = 1.959964 * Math::Sqrt(sigma2 / N);
richTextBox1->AppendText
("\n95%信頼区間= [ " + (p_hat - hfwidth95).ToString("F5")
+ ", " +
(p_hat + hfwidth95).ToString("F5") + " ]");
プログラムmleの実行形式のファイルmle.exeはファイル名mle.exeをクリックしてダウンロード・実行することができる。ただし、.NETFramework3.5対応である。プログラムを実行すると図１のフォームが表示される。
図１
全試行回数をＮに、Ｎ回の試行中、当該の結果（例えば、コインの表）の出た試行数をkに設定してCalcボタンを押すと計算が始まり、計算結果が図２のように表示される。
図２
プログラムのソースコードなどは圧縮形式のファイルmle.zipにまとめたので、ファイル名mle.zipをマウスの右ボタンでクリックすればダウンロードできる。ダウンロードしたファイルをマウスの右ボタンでクリックして表示されるメニュから「解凍」あるいは「展開」などを選んで解凍すればVisual Studioなどでプログラムを開くことができる。
次に、ベイズ的方法による分析のプログラムを説明する。最尤法においてフィッシャー情報行列を用いた区間推定を行うときは漸近定理に基づいているのでデータ数が十分に多いなどの条件がある。「十分に多い」ということの意味はデータに依存しており、100個で十分である場合があれば、300個でも不十分である場合もある。これに対して、ベイズ的方法はサンプル数が少ない場合でも理論的に適用可能である。ベイズ的方法については、「岡本安晴『統計学を学ぶための数学入門[下]』2009、培風館」の「9.12節 確率モデルによるデータ分析」における「ベイズ的分析」、「岡本安晴『データ分析のための統計学入門』2009、おうふう」の「6.4節 ベイズ的方法」、「岡本安晴『計量心理学』2006、培風館」の「付録H ベイズ的アプローチ」で説明している。
２項分布の場合、N回の試行において当該の事象がM試行生起したときの尤度関数は
で与えられるので、事前分布として一様分布をとれば事後分布は次式与えられる。
１変数の事後分布であるのでグリッド法により簡単に求めることができる。プログラムBayesBinでは事後分布を配列postPに以下のように求めている。
for (int i = 0; i <= 10000; i++)
sum
+= postP[i] = L(i/10000.0);
for (int i = 0; i <= 10000; i++){
postP[i]
/= sum;
if (maxV < postP[i]) maxV = postP[i];
}
ここで、L(p)は尤度関数である。事後分布から区間推定を行うが、これには確信区間（credible interval, central interval）と最高密度区間（Highest Density Region; HDR, Highest Posterior Density; HPD）がある。プログラムBayesBinでは、HPDをChen-Shao HPD estimation algorithmによって求めている。プログラムを以下に示す。
double
ckIntvl = 1.0;
int iL0 =
0, iR0 = 10000;
int ipos
= 0;
double
AlphaL = postP[ipos];
while
(AlphaL <= 0.05) {
int iposR = 10000;
double AlphaR = postP[iposR];
while (AlphaR < (0.05 - AlphaL)) AlphaR +=
postP[--iposR];
if (AlphaR > (0.05 - AlphaL)) iposR++;
if ((iposR - ipos) < (iR0 - iL0)){
iL0
= ipos; iR0 = iposR;
}
AlphaL
+= postP[++ipos];
}
int iL1 =
0, iR1 = 10000;
ipos = 10000;
double
AlphaR = postP[ipos];
while
(AlphaR <= 0.05) {
int iposL = 0;
double AlphaL = postP[iposL];
while (AlphaL < (0.05 - AlphaR)) AlphaL +=
postP[++iposL];
if (AlphaL > (0.05 - AlphaR)) iposL--;
if ((ipos - iposL) < (iR1 - iL1)){
iL1
= iposL; iR1 = ipos;
}
AlphaR
+= postP[--ipos];
}
double IntvlL = (iL0 + iL1) / 20000.0, IntvlR = (iR0 + iR1) / 20000.0;
プログラムBayesBinの実行形式のファイルBayesBin.exeはファイル名BayesBin.exeをクリックしてダウンロード・実行することができる。ただし、.NETFramework3.5対応である。プログラムを実行すると図3のフォームが表示される。
図３
全試行回数をＮに、Ｎ回の試行中、当該の結果（例えば、コインの表）の出た試行数をMに設定してDrawボタンを押すと計算が始まり、計算結果が図４のように表示される。
図４
ベイズ的分析において基礎となる結果は事後分布であるが、事後分布の代表値として中央値（Mode）と95％HPDが表示されている。NとMの値は画面左上のテキストボックス（図５）に設定し直して、Drawボタンをクリックすれば、設定した値で再計算が行われ計算結果が表示される。
図５
表示画面はSaveボタンをクリックすることによりファイルに保存することができる。
プログラムのソースコードなどは圧縮形式のファイルBayesBin.zipにまとめたので、ファイル名BayesBin.zipをマウスの右ボタンでクリックすればダウンロードできる。ダウンロードしたファイルをマウスの右ボタンでクリックして表示されるメニュから「解凍」あるいは「展開」などを選んで解凍すればVisual Studioなどでプログラムを開くことができる。