２項分布による簡単なベイズ分析
グリッド法
２項分布のようにパラメータが１個で値域が有界な場合は、ベイズ統計分析（ベイズ統計分析の基礎的説明のウェブサイトも用意している）における事後分布はグリッド法で簡単に求めることができる。
２項分布の場合、N回の試行において当該の事象が
試行生起したときの尤度関数は
で与えられる。したがって、事前分布として一様分布をとれば、事後分布は次式で与えられる。
尤度関数を次のスクリプト
def LF(p, N, k):
return (p ** k) * ((1.0 - p)
** (N - k))
で与え、グリッド法によって事後確率を得るためのクラス型PostPを以下のように宣言する。
n_grid = 100000
class PostP:
def __init__(self, N, k):
self.N = N
self.k = k
self.value = []
sum
= 0.0
for
i in range(0, n_grid + 1):
v = LF(i / n_grid, self.N, self.k)
self.value.append(v)
sum += v
for
i in range(0, n_grid + 1):
self.value[i] =
self.value[i] / sum
def d_p(self, i):
return self.value[i]
クラスPostPのメソッドd_pによって i 番目のグリッドの値が得られる。グリッドは、０から１までの区間 [0, 1] がn_grid = 100000 個に区切られているものである。０番目のグリッドの点が確率p の値０に、n_grid = 100000 番目のグリッドの点が確率p の値１に対応している。このクラス型を用いた例を以下に示す。このプログラム例は、実行すると分析結果が図１のように表示される。
図１
Ｎが総試行数、ｋが当該事象の生起数であり、成功確率 p のＭＡＰ（Maximum A Posteriori；事後分布のモード）推定値と９５％確信区間（95% CI）がグラフの上部に表示されている。Ｎとｋの値は、プログラムを実行すると図２のように聞いてくるので、計算したい値を設定すればよい。
図２
ｋの値を設定してEnterキーを押すと、図３のように計算結果が書き出され、図１のグラフが表示される。
図３
図１のフォームの右上のＸ印をクリックして閉じると、次のＮの値の設定が求められる（図４）。
図４
次の計算用のＮとｋの値を設定すれば、それに対応した計算が行われ、グラフが表示される。Ｎに負の値を設定すると、プログラムの実行終了となる。
プログラムのリストは以下に示すとおりである。なお、図１、図４では９５％確信区間（95% CI）が求められているが、９５％最高密度区間（95% HDR；Highest Density Region）（図５）を求めるプログラムも作成した。これは、このウェブサイトの後半に載せてある。これらのプログラムのスクリプトファイルは圧縮ファイルBinBysnPrg.zipにまとめた。
import
matplotlib.pyplot as plt
import numpy as np
def LF(p, N, k):
return (p ** k) * ((1.0 - p)
** (N - k))
n_grid = 100000
class PostP:
def __init__(self, N, k):
self.N = N
self.k = k
self.value = []
sum
= 0.0
for
i in range(0, n_grid + 1):
v = LF(i / n_grid, self.N, self.k)
self.value.append(v)
sum += v
for
i in range(0, n_grid + 1):
self.value[i] = self.value[i] / sum
def d_p(self, i):
return self.value[i]
while (True):
N = int(input('N = '))
if (N < 1):
break
k = int(input('k = '))
postP = PostP(N, k)
alpha = 0.05
L_p = -1
sum = 0.0
while (sum < alpha /
2.0):
L_p
+= 1
sum
+= postP.d_p(L_p)
#print('sum(L) = ', sum,
' L_p = ', L_p)
L_p = L_p / n_grid
R_p = n_grid + 1
sum = 0.0
while (sum < alpha /
2.0):
R_p
-= 1
sum
+= postP.d_p(R_p)
#print('sum(R) = ', sum,
' R_p = ', R_p)
R_p = R_p / n_grid
print('\nN = ', N, ' k = ', k)
print('MAP = ', k / N)
print('95% CI = [{0:.3f},
{1:.3f}]'.format(L_p, R_p))
plt.title(('MAP =
{0:.3f} for N = {1}, k =
{2}\n'.format(k / N, N, k)) +
('95% CI = [{0:.3f}, {1:.3f}]'.format(L_p, R_p)))
x = []
y = []
for i in range(0, n_grid +
1, n_grid // 1000):
x.append(i / n_grid)
y.append(postP.d_p(i))
plt.plot(x, y, 'b-')
plt.show()
９５％最高密度区間を求めるプログラム
以下のプログラムでは、図５に示すように９５％最高密度区間（95% HDR）が求められる。
図５
図５のグラフを表示するスクリプトは以下の通りである。
import
matplotlib.pyplot as plt
import numpy as np
def LF(p, N, k):
return (p ** k) * ((1.0 - p)
** (N - k))
n_grid = 100000
class PostP:
def __init__(self, N, k):
self.N = N
self.k = k
self.value = []
sum
= 0.0
for
i in range(0, n_grid + 1):
v = LF(i / n_grid, self.N, self.k)
self.value.append(v)
sum += v
for
i in range(0, n_grid + 1):
self.value[i] = self.value[i] / sum
def d_p(self, i):
return self.value[i]
while (True):
N = int(input('N = '))
if (N < 1):
break
k = int(input('k = '))
postP = PostP(N, k)
alpha = 0.05
temp_L = 0
temp_R = n_grid
L_p = temp_L
while True:
if
L_p == 0:
SumP = postP.d_p(L_p)
R_p = L_p
else:
SumP -= postP.d_p(L_p - 1)
while (SumP < 1 - alpha):
R_p += 1
if R_p > n_grid:
break
SumP += postP.d_p(R_p)
if R_p
> n_grid:
break
if
(temp_R - temp_L) > (R_p - L_p):
temp_L = L_p
temp_R = R_p
L_p
+= 1
if
L_p > n_grid:
break
L_p = temp_L / n_grid
R_p = temp_R / n_grid
print('\nN = ', N, ' k = ', k)
print('MAP = ', k / N)
print('95% HDR = [{0:.3f},
{1:.3f}]'.format(L_p, R_p))
plt.title(('MAP =
{0:.3f} for N = {1}, k =
{2}\n'.format(k / N, N, k)) +
('95% HDR = [{0:.3f}, {1:.3f}]'.format(L_p, R_p)))
x = []
y = []
for i in range(0, n_grid +
1, n_grid // 1000):
x.append(i / n_grid)
y.append(postP.d_p(i))
plt.plot(x, y, 'b-')
plt.show()