これはなに
- 二項分布で最尤推定/MAP推定/ベイズ推定をやります。
- 「イカサマコインの例で最尤推定とベイズ推定の違いを理解してみる」という記事にコメントしたものを多少編集したものです。
- 間違っている所があったらコメントください。
問題設定
- コインの表が出る確率を知りたい(0.7以上ならイカサマと判定する)
- コインを5回投げたら表3回・裏2回だった
- 表になる確率は0.75という事前知識がある
最尤推定/MAP推定/ベイズ推定の違い
| 手法 | 方法 | 事前確率の仮定 | パラメータ |
|---|---|---|---|
| 最尤推定 | 尤度関数を最大化するパラメータの値を決める | 仮定しない | 点推定 |
| MAP推定 | 尤度関数×事前分布を最大化するパラメータの値を決める | 仮定する | 点推定 |
| ベイズ推定 | パラメータの確率分布(事後分布)を決める | 仮定する | 事後分布 |
最尤推定について
1回のコイン投げをベルヌーイ分布で定式化します。
確率θで1(表)、確率1-θで0(裏)とすると、「5回試行して3回表が出る確率」は二項分布に従います。
これが尤度関数です。θを入力にとり、確率を出力します。
L(\theta) = {}_5C_3 \theta^3 (1-\theta)^2
θを0から1まで変化させた時の尤度関数は以下のように可視化できます。
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import beta, binom
sns.set()
theta = np.linspace(0, 1, 1000)
likelihood = binom.pmf(k=3, n=5, p=theta)
plt.plot(theta, likelihood)
plt.xlabel("Theta")
plt.ylabel("Likelihood")
plt.show()
尤度関数を最大化するθを求めます。
以下の式を解くと尤度関数を最大化するθがわかります。
\begin{eqnarray}
\frac{dL(\theta)}{d\theta} &=& 0 \\
10\theta^2(5\theta-3)(\theta-1) &=& 0 \\
\therefore \theta &=& \frac{3}{5}
\end{eqnarray}
尤度関数を最大化するθは対数尤度関数も最大にします。
\ln L(\theta) = \ln ({}_5C_3) + 3\ln (\theta) + 2\ln (1-\theta)
すなわち、以下の式を解いても良いです。
\begin{eqnarray}
\frac{d\ln L(\theta)}{d\theta} &=& 0 \\
\frac{3}{\theta}-\frac{2}{1-\theta} &=& 0 \\
\therefore \theta &=& \frac{3}{5}
\end{eqnarray}
MAP推定について
θの事前分布p(θ)を仮定します。尤度が二項分布に従う時、共役事前分布はベータ分布です。
ベータ分布は2つのパラメータa,bを持ちます。a,bは正の実数です。期待値はa/(a+b)です。
平均が0.75という情報が有るので、a/(a+b)=0.75とおくとa=3bになります。
というわけで、正の実数bに対してBeta(3b, b)という事前分布を仮定します。
b=1, 5, 10, 50の時のBeta(3b, b)の確率密度関数を可視化します。
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import beta, binom
sns.set()
bs = [1, 5, 10, 50]
x = np.linspace(0, 1, 1000)
for b in bs:
plt.plot(x, beta.pdf(x,3*b,b), label=f"b={b}")
plt.xlabel("x")
plt.ylabel("Probability density")
plt.legend()
plt.show()
事後確率を最大化するθを求めます。
ベイズ則から、事後確率=尤度x事前確率/周辺分布です。
周辺分布はθに依存しないので、事後確率を最大化するθは尤度x事前確率を最大化します。
L(\theta) \times Beta(\theta, 3b, b) = C\theta^{3+3b-1}(1-\theta)^{2+b-1} (Cは\thetaに依存しない定数)
同様に、事後確率を最大化するθはlog(尤度x事前確率)を最大化します。
以下の式を解くと事後確率を最大化するθがわかります。
\begin{eqnarray}
\frac{\partial \ln \{L(\theta) \times Beta(\theta, 3b, b)\}}{\partial \theta} &=& 0 \\
\frac{3+3b-1}{\theta}-\frac{2+b-1}{1-\theta} &=& 0 \\
\therefore \theta &=& \frac{3b+2}{4b+3}
\end{eqnarray}
b=[1, 100]の時のθの値を可視化します。
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import beta, binom
sns.set()
b = np.arange(1,100)
theta = (3*b+2)/(4*b+3)
plt.plot(b, theta)
plt.xlabel("b")
plt.ylabel("theta_map")
plt.show()
ベイズ推定について
θの事前分布p(θ)をBeta(3b, b)で仮定する所まではMAP推定と同じです。
今回は事前分布にBeta(6,2)を仮定します。
事後確率を最大化するθではなく、事後確率(θの関数)を直接計算します。
事後確率=尤度x事前確率/周辺分布の内、尤度と事前確率は計算済みです。
周辺分布p(D)を求めます。(p(5回試行して3回表が出る)の事です)
\begin{eqnarray}
p(D) &=& \int p(D, \theta)d\theta \\
&=& \int_{0}^{1} p(D|\theta)p(\theta)d\theta \\
\end{eqnarray}
尤度x事前分布p(D|θ)xp(θ)、周辺分布p(D)を計算し、事後分布p(θ|D)を求めます。
ここで、B(a, b)はベータ関数で、θに依存しない事を利用しています。また、5C3=10です。
\begin{eqnarray}
p(D|\theta) &=& {}_5C_3 \theta^3 (1-\theta)^2 \\
p(\theta) &=& \frac{\theta^{6-1} (1-\theta)^{2-1}}{B(6, 2)} \\
p(D|\theta)p(\theta) &=& \frac{10}{B(6, 2)} \theta^8 (1-\theta)^3 \\
&=& \frac{10}{B(6, 2)} (\theta^8 - 3\theta^9 + 3\theta^{10} - \theta^{11}) \\
p(D) &=& \int_{0}^{1} p(D|\theta)p(\theta)d\theta \\
&=& \frac{10}{B(6, 2)} \left[\frac{1}{9}\theta^{9} - \frac{3}{10}\theta^{10} + \frac{3}{11}\theta^{11} - \frac{1}{12}\theta^{12}\right]^1_0 \\
&=& \frac{1}{198 B(6, 2)} \\
p(\theta|D) &=& \frac{p(D|\theta)p(\theta)}{p(D)} \\
\therefore p(θ|D) &=& 1980\theta^{8}(1-\theta)^{3}
\end{eqnarray}
事前分布と事後分布を可視化します。
import numpy as np
import matplotlib.pyplot as plt
import seaborn as sns
from scipy.stats import beta, binom
sns.set()
theta = np.linspace(0, 1, 1000)
plt.plot(theta, beta.pdf(theta, 6, 2), label="Prior")
plt.plot(theta, 1980*theta**8*(1-theta)**3, label="Posterior") # 丸め誤差があるので本当は愚直に計算しないほうが良い
plt.xlabel("theta")
plt.ylabel("Probability density")
plt.legend()
plt.show()
事後分布をθで積分すると累積分布関数が求まります。
イカサマコインの閾値が0.7なので、θを[0.7, 1]で積分すればイカサマである確率を計算できそうです。
\begin{eqnarray}
\int^{1}_{0.7} p(\theta|D)d\theta &=& 1980\left[\frac{1}{9}\theta^{9} - \frac{3}{10}\theta^{10} + \frac{3}{11}\theta^{11} - \frac{1}{12}\theta^{12}\right]^{1}_{0.7} \\
&\sim& 0.507
\end{eqnarray}
約50.7%の確率でイカサマである事がわかりました。