コイン投げを例題にPythonを使ってベイズ統計学の勉強

概要

次のようなコイン投げの問題をベイズ統計学の観点から考えていきます。計算や結果の確認に Python を利用していきます。なお、ベイズの定理の基本については、理解していることを前提とします（ベイズの定理の基本については、「天気の子」を題材にゼロから学ぶ「ベイズの定理」を参照ください）。

あるコインについて「イカサマなコイン」かどうかを調べたい。このコインを3回投げたとき「裏」「表」「表」がでた。

コインの表がでる確率（スタンダードな統計学）

スタンダードな統計学（ネイマン・ピアソン統計学）では、次のように扱います。

あるコインについて「イカサマなコイン」かどうかを調べたい。このコインを3回投げたとき「裏」「表」「表」がでた。そのコインにおいて表がでる確率 $P(表)$ を求めよ。

表がでる確率 $P(表)$ は、次のように求めることができます。

$$P(表)=\frac{2}{3}$$

ところで、この $P(表)$ は、6回投げて「裏・表・表・裏・表・表」となった場合でも同じく $2/3$ となります。当然、3回よりも6回投げたときに $P(表)=2/3$ となったほうが確信をもって「イカサマなコイン」といえるのですが、このようなことを考慮してイカサマ判定するためには仮説検定という手法が必要となります。仮説検定は、ここでの本題ではないので省略します。

コインの表がでる確率の確率分布（ベイズ統計学）

ベイズ統計学では、母数に関する確率分布（＝母数を確率変数とした確率密度関数）を扱います。母数とは、母集団の性質を決めるパラメータのことで、このコインの例であれば「表がでる確率」、正規分布の推定であれば「平均 $\mu$」や「標準偏差 $\sigma$」となります。

コイン投げの問題を、ベイズ統計学では次のように扱います。

あるコインについてイカサマなコインかどうかを調べたい。このコインを3回投げたとき「裏」「表」「表」がでた。このコインについて「表がでる確率 $\theta$」の確率分布を求めよ（表が出る確率 $\theta$ を確率変数としたときの確率密度関数 $f(\theta)$ を求めよ）。

「表がでる確率 $\theta$」ではなく「表がでる確率 $\theta$ の確率分布（確率密度関数 $f(\theta)$ ）」を求めるところがベイズ的なポイントとなります。「表がでる確率 $\theta$」は $0.0$ ～ $1.0$ の範囲をとるので、求めるべき確率密度関数 $f(\theta)$ は、次のようなグラフ上に何らかの形として描かれるはずです（無論、現時点では、どのようになるかは分かりません）。

ところで、この確率密度関数 $f(\theta)$ は、$\theta$ がとりうる範囲 $[0.0, 1.0]$ で積分したときに $1.0$ となる性質（面積が $1.0$ となる性質）を持ちます。

$$ \int_{0.0}^{1.0} f(\theta)\, d\theta = 1.0 $$

話を戻します。ベイズ統計学において、確率密度関数 $f(\theta)$ は、ベイズの定理を使って求めていきます。ベイズの定理とは、次のようなものです。

$$P(\theta|y) = \frac {P(y|\theta)\, P(\theta)}{P(y)} $$

• $P(\theta | y)$ は「事象 $y$（表or裏）が観測されたという結果」があるときに、その原因として推定される $\theta$ の確率分布で、これを事後分布といいます。
• $P(\theta)$ は、事象が観測される前段階において推定されている $\theta$ の確率分布で、これを事前分布といいます。事前分布に対する適切な情報がない場合、主観により任意に（恣意的に）設定することができます。
• $P(y|\theta)$ は、表がでる確率を $\theta$ と仮定したときに、事象 $y$ が観測される確率になります。今回の場合、観測された事象 $y$ が「表」ならば $P(表|\theta)=\theta$ となり、「裏」ならば $P(裏|\theta)=(1-\theta)$ となります。尤度（ゆうど）ともいいます。
• $P(y)$ は、事前分布のもとで事象 $y$ が観測される全確率になります。

以上に基づいて、コインの表がでる確率 $\theta$ の確率分布 $f(\theta)$ を求めていきます。

準備：事前分布の設定

はじめに、「裏・表・表という観測結果が得られる」ということは知らない前提で事前分布 $P(\theta)$ を決めていきます。いまは「一切の先入観なしでイカサマなコインであるかを調べたい」という状況から事前分布 $P(\theta)$（表がでる確率 $\theta$ の確率密度関数）は $1$ であるとします（理由不十分の原則）。つまり、表がでる確率 $\theta$ が、$\theta=0.5$ である可能性も、$\theta=0.0$ である可能性も、$\theta=1.0$ である可能性もすべてが等しい、というところをスタートにします。

実は、ここがベイズ統計学の特徴的なところでもあるのですが、別に事前分布を $1$ で与える必要はまったくなく、コインの外観と手に持った感触から「表と裏は均等にでそうだ」と思えば、$\theta$ が $0.5$ にピークを持つような事前分布 $P(\theta)$ を与えることもできます。このケースについては後ほど扱います。

とりあえず、ここでは事前分布 $P(\theta)=1$ とします。この事前分布は、次のように描けます（このグラフを描くためのPythonコードは最後に掲載しています）。$f(\theta)$ を積分すると（面積を計算すると）$1.0$ になり、確率密度関数として適切なものであることが確認できると思います。

1回目に「裏」を観測

1回目のにコインを投げたとき「裏」という事象が観測されました。この情報を使って、1回目の試行後に推定される $\theta$ の確率分布（＝事後分布） $P(\theta | y)$ を次の「ベイズの定理」により計算していきます。

$$P(\theta|y) = \frac {P(y|\theta)\, P(\theta)}{P(y)} $$

前節で決めた事前分布 $P(\theta)$ は $1$、また $P(裏|\theta)$ は、表がでる確率を $\theta$ としたときの「裏のでる確率」なので $(1-\theta)$ になります。また、分母の $P(裏)$ は確率密度関数 $f(\theta)=1$ のときに裏がでる確率で、次のように簡単な積分を使って求められます。

$$ \int_{0.0}^{1.0}(1-\theta)\cdot f(\theta)\,d\theta = \int_{0.0}^{1.0}(1-\theta)\cdot 1\,d\theta = \frac{1}{2} $$

積分計算には SymPy が利用できます（実行すると 1/2 が得られます）。

定積分

import sympy
t = sympy.Symbol('t')
sympy.integrate((1-t)*1, (t, 0, 1))

これをベイズの定理の右辺に代入すると、事後分布 $P(\theta | y)$ は次のようになります。

\begin{align}
P(\theta|y) & = \frac {P(y|\theta)\, P(\theta)}{P(y)} \\
            & = \frac {(1-\theta )\cdot 1 }{\frac{1}{2}} \\
            & = 2-2\theta 
\end{align}

事後分布 $P(\theta | y)=2-2 \theta$ を図示すると次のようになります。

試行の結果、「裏」という事象を観測したため、表が出る確率 $\theta$ が 低くなるような分布に変化していることがわかると思います。この分布に従えば、$\theta$ が $0.0$～$0.1$ にある確率は $19\%$、 $0.9$～$1.0$ にある確率は $1\%$ と計算できます。

$$\int_{0.0}^{0.1} f(\theta)\,d\theta = \int_{0.0}^{0.1}(2-2\theta)\,d\theta = 0.19 $$
$$ \int_{0.9}^{1.0} f(\theta)\,d\theta = \int_{0.9}^{1.0}(2-2\theta)\,d\theta = 0.01$$

import sympy
t = sympy.Symbol('t')
print( f'p1 = {sympy.integrate( 2-2*t, (t, 0.0, 0.1)):>5.1%}' )
print( f'p2 = {sympy.integrate( 2-2*t, (t, 0.9, 1.0)):>5.1%}' )

実行結果

p1 = 19.0%
p2 =  1.0%

2回目に「表」を観測

2回目のコイン投げをしたら今度は「表」という事象が観測されました。1回目と同様にベイズの定理で事後分布を求めていきますが、この際、事前分布 $P(\theta)$ には、1回目の事後分布 $(2-2 \theta)$ を利用します。つまり、1回目の観測により得られた情報を利用して、2回目の試行後の確率分布 $P(\theta|y)$ を推定するわけです。これをベイズ更新といいます。

$$P(\theta|y) = \frac {P(y|\theta)\, P(\theta)}{P(y)} $$

ベイズ更新により $P(\theta)$ は $(2-2 \theta)$です。$P(表|\theta)$ は、表がでる確率を $\theta$ としたときの「表のでる確率」なので $\theta$ になります。また、$P(表)$ は確率密度関数 $f(\theta)=(2-2\theta)$ のときに表がでる確率で、次のような計算で求まります。

$$ \int_{0.0}^{1.0}\theta\cdot f(\theta)\,d\theta = \int_{0.0}^{1.0}\theta\cdot(2-2\theta)\,d\theta = \frac{1}{3} $$

定積分

import sympy
t = sympy.Symbol('t')
sympy.integrate(t*(2-2*t), (t, 0, 1)) # 1/3

これらをベイズの定理の右辺に代入すると、2回目の試行後に推定される $\theta$ の確率分布、つまり、「裏」「表」が観測されたときの $\theta$ の確率分布が計算できます。

\begin{align}
P(\theta|y) & = \frac {P(y|\theta)\, P(\theta)}{P(y)} \\
            & = \frac {\theta \cdot (2-2\theta) }{\frac{1}{3}} \\
            & = 6\theta -  6\theta^2  
\end{align}

事後分布 $P(\theta | y)=6\theta - 6\theta^2$ を図示すると次のようになります。

確率密度関数 $f(\theta)$ が、$\theta=0.5$ でピークを持つ形に変化したことが分かります。この確率密度関数から言えることは、表がでる確率が 50% 付近である可能性が高いということです。これは、私たちの直感的な理解とも一致します（コインを2回投げて表と裏が1回ずつでたら、表と裏がでる確率はおそらく同じだろう・・・と思いますよね）。

3回目に「表」を観測

つづいて、3回目のコイン投げをしたら今度も「表」という事象が観測されました。
例によって、ベイズの定理から3回目の試行後に推定される $\theta$ の確率分布（＝事後分布）を求めるわけですが、そこで利用する事前分布 $P(\theta)$ は、2回目の事後分布 $6\theta-6\theta^2$ です（ベイズ更新）。

$P(表|\theta)$ は、表がでる確率を $\theta$ としたときの「表のでる確率」なので $\theta$ になります。また、$P(表)$ は確率密度関数 $f(\theta)=6\theta-6\theta^2$ のときに表がでる確率で、次のように積分により求められます。

$$ \int_{0.0}^{1.0}\theta\cdot f(\theta)\,d\theta = \int_{0.0}^{1.0}\theta\cdot(6\theta-6\theta^2)\,d\theta = \frac{1}{2} $$

定積分

import sympy
t = sympy.Symbol('t')
sympy.integrate(t*(6*t-6*t**2), (t, 0, 1)) # 1/2

これらをベイズの定理の右辺に代入すると、3回目の試行後に推定される $\theta$ の確率分布、つまり「裏」「表」「表」が観測されたときの $\theta$ の確率分布が計算できます。

\begin{align}
P(\theta|y) & = \frac {P(y|\theta)\, P(\theta)}{P(y)} \\
            & = \frac {\theta \cdot (6\theta-6\theta^2) }{\frac{1}{2}} \\
            & = 12\theta^2 -  12\theta^3 
\end{align}

事後分布 $P(\theta | y)=12\theta^2 - 12\theta^3$ を図示すると次のようになります。

この確率密度関数から分かるのは、表がでる確率 $\theta$ は $2/3$ である可能性にピークを持つような分布となっているということです。コインを3回投げて「裏」「表」「表」が観測されたことから考えられる表の出現確率に関する直感的な理解,スタンダードな統計学との結果とも一致します。

表と裏がでた順番は最終的な確率分布に影響するのか？

ここまで、コインを投げたときに「裏」「表」「表」という順で結果がでた事例を扱いました。この順番が変わったとき、最終的な事後分布は影響をうけるのでしょうか？

答えは「影響されない」です。たとえば「表」「表」「裏」という順であっても最後に求まる事後分布は、先ほどと同じ $12\theta^2 - 12\theta^3 $ となります。これは、数式から確認することができます。3回目の事後分布を計算する際の事前分布には「2回目の事後分布」を使いました。2回目の事後分布を計算する際の事前分布には「1回目の事後分布」を使いました。

つまり、3回目の事後分布は次のように計算できることになります。

$$ P(\theta|y) = f(\theta) = \frac {\theta \cdot \theta \cdot (1-\theta) \cdot 1}{Z}=\frac{\theta^2-\theta^3}{Z} $$

式の分子の最初にでてくる $\theta$ は3投目で「表」がでたこと、次の $\theta$ は2投目で「表」がでたこと、次の $(\theta-1)$ は1投目で「裏」がでたこと、次の $1$ は観測前の事前確率として $1$ を設定したことに対応します。1投目から3投目までの結果に対応する部分は、掛け算なので順番が入れ替わっても結果に影響しません（これを逐次合理性といいます）。

ところで、$Z$ は $f(\theta)$ が適切な確率密度関数（つまり、$0.0$ から $1.0$ の範囲で積分したときに $1.0$ ）となるように決める定数となります。つまり、次式を満たすような定数 $Z$ です。

$$ \int_{0.0}^{1.0}\frac{\theta^2-\theta^3}{Z}\,d\theta = 1.0$$

SymPyで $Z$ を計算してみます。

import sympy
t = sympy.Symbol('t')
Z = sympy.Symbol('Z')
sympy.solve(sympy.integrate((t**2-t**3)/Z,(t,0,1))-1, Z) # 1/12

$Z=1/12$ と求まりました。これを先ほどの式に代入すれば、前節と同じ結果が得られます。つまり、ちまちまと計算しなくとも一気に計算ができるということが分かります（積分計算も1回だけでOKなのです）。

$$ P(\theta|y) = f(\theta) = \frac{\theta^2-\theta^3}{Z} = \frac{\theta^2-\theta^3}{\frac{1}{12}} = 12\theta^2-12\theta^3$$

一般化できないのか？

できます。コインを投げたとき、表が $a$ 回、裏が $b$ 回観測されたときの「表がでる確率 $\theta$ 」の確率分布（確率密度関数）は、事前分布を $P_{0}(\theta)$ とすれば次のようになります。

$$ P(\theta|y) = \frac{\theta^{\, a}\cdot(1-\theta)^b\cdot P_{0}(\theta)}{Z} $$

ただし、$Z$ は、次式を満たす定数です。

$$ \int_{0.0}^{1.0} \frac{\theta^{\, a}\cdot(1-\theta)^b\cdot P_{0}(\theta)}{Z}\,d\theta = 1.0$$

補足：ベータ分布・ベータ関数

じつは、上記の $P(\theta|y)$ はベータ分布という確率密度関数の形となっています。ベータ分布は、次のように定義される確率密度関数です。

$$ f(x;\alpha,\beta) = \frac{x^{\alpha-1}\cdot (1-x)^{\beta-1}}{B(\alpha,\beta)} $$

ここで、$B(\alpha,\beta)$ はベータ関数というものになります。このベータ関数の値は scipy.special.beta でサラリと求めることができます。

表が4回・裏が2回でたときの θ の確率分布

一般化した式を利用して「表が4回・裏が2回」でたときの $\theta$ の確率分布を求めます。観測前の事前分布は $1$ とします。

$$ P(\theta|y) = \frac{\theta^{\, a}\cdot(1-\theta)^b\cdot P_{0}(\theta)}{Z} = \frac{\theta^{\, 4}\cdot(1-\theta)^2\cdot 1}{Z} $$

$P(\theta|y)$ が確率密度関数の要件を満たすように $Z$ を求めます。

import sympy
t = sympy.Symbol('t')
Z = sympy.Symbol('Z')
sympy.solve(sympy.integrate((((t**4)*(1-t)**2))/Z,(t,0,1))-1, Z) # 1/105

求められた $Z=1/105$ を式に代入すれば、表がでる確率 $\theta$ の確率分布は次のようになります。

$$ P(\theta|y) = 105\cdot\theta^{\, 4}\cdot(1-\theta)^2 $$

なお、$Z$ はベータ関数から次のように求めることもできます。

import scipy
scipy.special.beta(4+1, 2+1) # 0.0095238... = 1/105

確率分布を図示すると次のようになります。

「3投中2回の表」のときと同じく確率密度のピークは $\theta=2/3$ にあります。ただし、先に示した「3投中2回の表」のときのグラフと比較すると「ばらつき」が小さくなっていることが分かると思います。例えば、$\theta$ が $0.0$ から $0.1$ となる確率はかなり低くなっています。

スタンダードな統計学では「3投中2回の表」でも「6投中4回の表」でも、結果として得られる確率は $2/3$ で違いはありませんが、ベイズ統計学では確率分布の違いとして表れてきます。なんらかの意思決定の資料と考えた場合、ベイズのほうが有用ですよね。

事前分布を 1 以外に設定できるのか？

ここまで、観測前の事前分布は $1$ としてきましたが、それ以外ものを設定することもできます。ただし、それは積分して $1.0$ となる確率密度関数でないと計算がおかしくなります。

例えば、コインの外観と手に持った感触から「表と裏は均等に出現しそうだ」と思えば、観測前の事前分布 $P(\theta)$ として $6\theta - 6\theta^2$ を設定することもできます。これは、すでに示していますが、以下のような確率分布になります。

図から分かるように「表と裏は均等に出現しそうだ・・・」ということを仮定に推定をスタートすることもできるのです。

ところで、左右対称で中央にピークを持つ分布は、正規分布など色々とあるなかで $6\theta-6\theta^2$ を使ったことには理由があります。$6\theta-6\theta^2$ は $\alpha=2$、$\beta=2$ のベータ分布なのです。

$$ 6\theta-6\theta^2 = \frac{\theta^{\,2-1}\cdot(1-\theta)^{2-1}}{1/6} = f(\theta\,;\alpha=2,\beta=2)$$

このコインの例の場合、事前確率をベータ分布の形にしておくと、そのあとの事後分布もベータ分布の形になるので、$Z$ を求めることが簡単になります（Zの計算にベータ関数が利用できます）。

なお、$1$ も $\alpha=1$、$\beta=1$ のベータ分布といえます。

$$ f(\theta\,;\alpha=1,\beta=1) = \frac{\theta^{\,1-1}\cdot(1-\theta)^{1-1}}{B(1,1)} = \frac{1\cdot 1}{1}=1 $$

確率密度関数のグラフを描くコード

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
from IPython.display import display

fp= FontProperties(fname='ipaexg.ttf', size=11);

xRange=(0,1)  # グラフ描画範囲 x
yRange=(0,2)  # グラフ描画範囲 y

plt.figure(figsize=(5,3), dpi=120)

f = lambda t : 6*t - 6*t**2 # ここを書き換え
x = np.linspace(xRange[0], xRange[1], 100)
y = f(x)

plt.plot(x, y, color='orange')

plt.xlim(xRange[0], xRange[1])
plt.xticks(np.linspace(xRange[0], xRange[1], 11))  
plt.xlabel( r'確率変数 $\theta$（表がでる確率）' , fontproperties=fp )

plt.ylim(yRange[0], yRange[1])
plt.yticks(np.linspace(yRange[0], yRange[1], 5))  
plt.ylabel(r'確率密度 $f\,(\theta)$', fontproperties=fp )

plt.savefig('test.svg',format = 'svg'); # SVG形式で出力

plt.show()

参考資料

• [1] 涌井貞美 著「図解・ベイズ統計「超」入門」おすすめ
• [2] 一石賢 著「意味がわかるベイズ統計学」
• [3] 小島寛之 著「完全独習ベイズ統計学入門」