ベイズ推定 Python 連携

ベイズ推定の概要

ベイズ推定とは、ベイズ定理の考え方に基づく統計的推定のことです。事前の確率分布と観測値に基づいて、事後の確率分布を求めます。観測値が多ければ多いほど、正確な結果に導く可能性が高まるため、データ処理性能が高まる近年、改めて注目を浴びる理論です。

事後の確率分布を事後分布と呼び、以下の計算式で求めることができます。

p(θ|x)=(p(x|θ) × p(θ)) ÷ (p(x))

事後分析=(尤度関数 × 事前分布) ÷ 規格化定数

ざっと用語を説明すると以下の通りですが、細かくは以後の流れを見ながら確認しましょう。

要素	式	説明
事後分布	p(θ|x)	事象xが発生した状況下で事象θが起こる確率
尤度関数	p(x|θ)	事象θが発生した状況下で事象xが起こる確率
事前分布	p(θ)	事象θが発生する確率
規格化定数	p(x)	事象xが発生する確率

ケース

100人の被験者を50人ずつの2グループに分け、それぞれに本物の薬と偽薬(プラセボ)を投与し、薬の効果を確かめる試験を実施しました。結果、本物の薬を投与した50人のうち40人に効果が確認できました。一方のプラセボを投与したグループでは、50人のうち10人に効果が確認できました。

この薬の効果をベイズ推定の考え方を用いて確認します。「事前分布」→「尤度」→「事後分布(規格化前)」→「事後分布」の順に値を求めます。

本物の新薬を投与したグループに所属する確率(p)と、偽薬を投与したグループに所属する確率(1−p)は、それぞれ0.5です。また、薬の効果が確認できた場合、その被験者が本物の薬を投与したグループに所属する確率は、0.1刻みで0.1～0.9の離散値を取ることにします。

事前分布

観測結果が得られるまで、実際の分布を描くことはできません。そのため、全ての事前確率は同じであることを前提に、以下のような一様分布を決定します。

Yellowfinの棒グラフを使って表すと以下の感じです。

尤度

ある前提条件に従って結果が出現する場合に、逆に観察結果からみて前提条件が「何々であった」と推測する尤もらしさ (出所：Wikipedia)を示すのが尤度です。
例えば成功率0.9の尤度は、50回の投薬のうち、40回成功する(効果が出る)という条件付き確率を意味し、以下の式で求めることができます。

p(k=40|n=50,p=0.9)=\begin{pmatrix} 50\\40 \end{pmatrix}(0.5)^{40}(1-0.5)^{10}≒0.01518

同様に、P=0.1から0.8の値も以下のPythonのコードで求めます。算出した値は、後続処理のためにdata変数に配列として取り込みます。

尤度関数

from scipy.stats import binom
data = []
for i in range(1,10):
    p = binom.pmf(40, 50, i/10)
    data.append(p)
print(data)

#実行結果
[3.5817219285889007e-31,1.212736553318045e-19,3.5277464245644363e-13,
7.508945087955164e-09,9.123615791750718e-06,0.001439848174336955,
0.038618990683862445,0.1398190051743154,0.015183334117262373]

Yellowfinの棒グラフでは以下のような感じになりました。

事後分布(規格化前)

尤度と事前分布を掛け合わせて足したものが規格化前の事後分布です。
規格化前の事後分布の数値の総和が規格化定数で、0.02146です。

事後分布は発生確率を表すものです。次のステップで、事後分布(規格化前)の各数値を規格化定数で割って、総和が1になるようにします。

事後分布

事後分布(規格化前)の各数値を規格化定数で割って、事後分布の数値を求めます。

p=0.1から0.9までの数値求めるPythonコードは以下になります。

事後分布

for i in range(0,9):
    print(data[i]/sum(data))

#実行結果 ↓
1.836118444628096e-30
6.216920236739975e-19
1.8084486755970698e-12
3.849353146497846e-0#8
4.67709095539237e-05
0.007381175432023921
0.19797472422851276
0.7167621033362965
0.07783518759827292

このデータをYellowfinに取り込んで可視化します。

P>0.5である確率は99.9%以上で、薬の効果は非常に高いことが分かります。

新薬を投与したグループの結果で、最も高い数値を示すのがp=0.8の0.7168でしたこの結果から、薬の効果が確認できた被験者が本物の薬を投与したグループの所属する確率が、プラセボを投与したグループに所属する確率の4倍である可能性が最も高いことが分かります。新薬の開発でこのような効果が確認できれば、直ぐにでも厚労省の承認が下りそうですね。知らんけど。

ダッシュボードの作成

作成したチャートをダッシュボードに配置してみます。

最後に

ベイズ推定は、事前の確率分布と観測値に基づいて、事後の確率分布を求めるため、観測値が多ければ多いほど、正確な結果に導く可能性が高まることを冒頭で説明しました。例えば、被験者の数を単純に10倍した数値で事後分布を求めると、以下のようにp=0.8の値がより突出することが確認できます。観測結果が大きければ大きいほど、正確な数値に近づく可能性が高まることを意味します。サイズの大きなデータを処理できる昨今、ベイズ推定(定理)の重要性が見直されていることが、この結果からも読み取ることができます。

今回は簡単なモデルを使って試してみました。今後、実データを使った検証も行ってみたいと思っています。

では皆様、良いデータ分析を! Cheers, geeks!!

参考情報