復習がてら、ベイズ統計の解説とPythonでの実装を行います。今のところ以下のように3つに分けて書いていこうと思います。
第1回：ベイズ推定とは
第2回：共役事前分布について
第3回：MCMC（Markov Chain Monte Carlo）による事後分布のサンプリング
今回は第1回ということで「ベイズ推定とは何か」についてです。
前提知識と記法
ベイズ統計“再入門”なので多少の統計の知識を仮定します（分布とか）。 高校数学でやる基本的な確率や微分の計算は既知とします。 Pythonはちょっと触ったことのある程度で大丈夫です。記法など
$X$と$Y$を事象とします。$P(Y)\neq 0$のとき $$ P(X\mid Y)=\frac{P(X\cap Y)}{P(Y)} $$で定義される確率を事象$Y$が起こった下で事象$X$が起きる条件付き確率と呼びます。ちなみに英語だと "probability of $X$ under the condition $Y$" なので、$X$と$Y$の登場する順番が同じで分かりやすい気がします。ともあれ、条件は後ろ。全事象$U$が有限で根源事象は同様に確からしいとし、事象$A$が起きる場合の数を$n(A)$で表すことにすると、条件付き確率は
$$
P(X\mid Y)=\dfrac{n(X\cap Y)\div n(U)}{n(Y)\div n(U)} = \frac{n(X\cap Y)}{n(Y)}
$$で計算できます（高校では有限個の事象しか扱わないのでこっちの式を条件付き確率の定義としていますね）。
コインを投げた結果から表の出る確率を推定する
この記事では次のコインの問題を考えます：裏表の出る確率が分からないコインを3回投げたとき、裏、表、表と出た。このコインの表が出る確率は何か？
この問いに対しての真の答えは神様にしか分かりません（たぶん）。しかし、統計学に基づいた予想（すなわち推定）をすることは可能です。
ところで、推定では「何を仮定するか」が重要になってきます。
極端な話「裏表の出る事象が同様に確からしい」と仮定すれば、表が出る確率は$\frac{1}{2}$です。これも推定といえば推定ですが、与えられたデータを全く利用していないので微妙です。
また、データから単純に表が出る確率を$\frac{2}{3}$だろうと推定する人もいると思います。結果だけ見ると実はこの推定は後述の「最尤推定」の結果と一致するのですが、この推定の問題点はデータ数が少ないと精度がわるくなる可能性があります。例えば仮に3回投げて3回とも表だった場合、同じ手法でこのコインの表が出る確率を推定すると1になります。「このコインは表しか出ない」を信じろと言われても結構無理があると思います。
他にも色々な推定が考えられていますが、この記事では「ベイズ推定」を紹介したいと思います。
最尤推定
まずはベイズ推定と比較するために最尤推定について解説します。$N$個のデータ$x_1,\ldots, x_N$が観測されたとします。このとき、パラメータ$\theta$を与えたときに、$N$個のデータが起こる条件付き確率
P(x_1,\ldots,x_N \mid\theta)
を求め、この値を最大にする$\theta=\hat{\theta}$こそが求めていたパラメータであろうと推定する手法が最尤推定です。
\hat{\theta} = \underset{\theta}{\operatorname{argmax}} P(x_1,\ldots,x_N \mid\theta)
最尤推定に基づいてコインの問題を考えます。コインの表が出る確率が$p$だったときにコインが裏、表、表と出る確率は
\begin{align}
P(T,H,H \mid p) &= (1-p)p^2\\
&= p^2 -p^3
\end{align}
です（Headは表、Tailは裏）。この値を$L(p)=p^2 -p^3$とおき、$L$を最大にする$p=\hat{p}$を求めます。これはふつうに微分で求めれられます。
$$
L'(p) = 2p - 3p^2 = p(2-3p)
$$よって、$L'(p)=0$となるのは$p=0, \frac{2}{3}$ですが、$L$のグラフを考えると$p=\frac{2}{3}$のときに$L$は最大値をとります。したがって、$\hat{p}=\frac{2}{3}$と推定できました。結果だけみると「コインが裏、表、表と出たから$\frac{2}{3}$」と考える方法と同じですね。
次にベイズ推定でコインの問題を考えます。
ベイズ推定
ベイズの定理と証明
ベイズ推定では次の**ベイズの定理**を利用します。$X, Y$を事象とし、$P(Y)\neq0$とする。このとき
$$
P(X\mid Y)=\dfrac{P(Y\mid X)P(X)}{P(Y)}
$$
証明
条件付き確率の定義より
\begin{align}
P(X\mid Y)P(Y) &= P(X\cap Y), \\
P(Y\mid X)P(X) &= P(Y\cap X)
\end{align}
なので$P(X\mid Y)P(Y)=P(Y\mid X)P(X)$が成り立つ。両辺を$P(Y)\neq0$で割れば求める式が得られる。■
ベイズ推定でコインの問題を考える
ベイズ推定では$N$個のデータ$x_1,\ldots, x_N$が観測されたときに、パラメータ$\theta$が起こる条件付き確率$$
P(\theta \mid x_1,\ldots,x_N)
$$を考えます（事後分布といいます）。ベイズの定理によって
$$
P(\theta \mid x_1,\ldots,x_N) = \frac{P(x_1,\ldots,x_N \mid\theta)P(\theta)}{P(x_1,\ldots,x_N)}
$$となるので、この確率を最大にするパラメータ$\theta=\hat{\theta}$こそが求めていた値であろうと推定する手法がベイズ推定（MAP推定）です。
$$
\hat{\theta} = \underset{\theta}{\operatorname{argmax}} \frac{P(x_1,\ldots,x_N \mid\theta)P(\theta)}{P(x_1,\ldots,x_N)}
$$ここで、$P(x_1,\ldots,x_N \mid\theta)$を尤度関数，$P(\theta)$を事前分布、$P(x_1,\ldots,x_N)$をエビデンスと呼びます。
この推定を用いてコインの問題を考えます。コインが裏、表、表と出たときに、コインの表が出る条件付き確率は
\begin{align}
P(p \mid T,H,H)
&= \frac{P(T,H,H \mid p)P(p)}{P(T,H,H)}\\
&= \frac{(p^2-p^3)P(p)}{P(T,H,H)}
\end{align}
となります（最尤推定のときに計算した結果を使いました）。分母の$P(T,H,H)$は計算できなくもないですが、ただの定数なので$P(p \mid T,H,H)$最大にする$\theta$を求めるのに影響しないので無視します。
分子の$P(p)$は「コインを投げて表がでる確率が$p$になる確率」を表します。つまり、この世のコインの表が出る確率を調べたときに得られる分布のことです。これを正確に知ることは難しいですね。
先ほども言いましたが、推定では「何を仮定するか」が重要です（ここが腕の見せ所）。おそらくですが、この世のコインの表が出る確率の平均は$0.5$で、標準偏差は$0.05$くらいではないでしょうか。そこで$P(p)$は平均$0.5$、標準偏差が$0.05$となるベータ分布に従うと仮定しましょう。
$$
P(p)=\frac{(1-p)^{49.5-1}p^{49.5-1}}{B(49.5,49.5)}=C(1-p)^{48.5}p^{48.5}
$$ここで$B(49.5,49.5)$は定数なので$C=\frac{1}{B(49.5,49.5)}$とおきました（なぜ正規分布ではなくてベータ分布を選んだかについては第2回で書こうと思います）。よって
\begin{align}
P(p \mid T,H,H)
&= \frac{(1-p)p^2 P(p)}{P(T,H,H)}\\
&= \frac{(1-p)p^2 C(1-p)^{48.5}p^{48.5}}{P(T,H,H)}\\
&= C'(1-p)^{49.5}p^{50.5}
\end{align}
となります（定数はまとめて$C'$と置きました）。$L(p)=(1-p)^{49.5}p^{50.5}$を最大にする$p$はこのまま普通に$L'(p)$を計算してもできますが
$$
\log L(p) = 49.5\log(1-p) + 50.5\log p
$$を微分して0になる値を求める方法でやってみましょう（尤度方程式）。
\begin{align}
\frac{d}{dp}\log L(p)=0 &\iff -\frac{49.5}{1-p} + \frac{50.5}{p} = 0\\
&\iff p = \frac{50.5}{100}
\end{align}
よって、$\hat{p}=\frac{50.5}{100} = 0.505$が求める推定値になります。若干表が出やすいという結果になりましたね。最尤推定よりも信頼できる結果になったと思います。
Pythonで計算
簡単ではありますが、今の計算をPythonで計算してみます。import numpy as np
from scipy.stats import beta
import matplotlib.pyplot as plt
%matplotlib inline
a,b = 49.5, 49.5
print(f"平均:{beta.mean(a, b)}, 標準偏差:{beta.std(a, b)}")
平均:0.5, 標準偏差:0.05
h = 2 #表の回数
t = 1 #裏の回数
def p_post(p,heads,tails):
return (p**heads) * ((1-p)**tails) * beta.pdf(p, a, b)
n = 10000
p = np.linspace(0,1, n)
p_map = np.argmax(p_post(p,h,t))/n
print(f"MAP推定値: {p_map}")
fig, ax = plt.subplots(1, 1)
plt.xlim(0.38, 0.62)
plt.plot(p, p_post(p,h,t))
plt.axvline(x=p_map, ymin=0, ymax=p_post(p_map,h,t), color="r")
MAP推定値: 0.5049
確かに0.505付近で最大値をとっていますね。
まとめと次回予告
今回はベイズ推定のなかでもMAP推定と呼ばれる推定について簡単に解説しました。ポイントは- 最尤推定は表が出る確率を仮定する
- ベイズ推定は事前分布を仮定する
の2つかなと思います。
次回は事前分布の選び方の解説と、もう少し複雑な例で実装をしたいと思います。