ベイズ推論イントロ

Pythonで体験するベイズ推論

Pythonで体験するベイズ推論 PyMCによるMCMC入門

Pythonで体験するベイズ推論 PyMCによるMCMC入門

  • 作者: キャメロンデビッドソン=ピロン,玉木徹
  • 出版社/メーカー: 森北出版
  • 発売日: 2017/04/06
  • メディア: 単行本(ソフトカバー)
  • この商品を含むブログを見る
から, ベイズ推論についての紹介

世界に司書と農家しか職業選択が無いときに司書と農家の人口比は1:20である.
Steveが内向的でおとなしい性格だという情報を得たときに, Steveは司書である確率はどのくらいなのであるのだろうか?

 P(A)をSteveが司書である確率とする.
このとき求めたいのはSteveが内向的な性格であるとわかった上で, Steveが司書であるという確率 P(A|X)である.
ベイズの定理より
 P(A|X) = \frac{P(X|A)P(A)}{P(X)}
 P(X|A)は司書であるときに内向的な性格であるという確率である. これは0.95とする.
 P(A)は司書である確率なので \frac{1}{21}である.
 P(X) = P(X かつ A) + P(X かつ not A) = P(X|A)P(A)+P(X|not A)P(not A)より,
 P(X|not A)は農家であるときに内向的な性格である確率で0.5とする.
すると P(X) = 0.95*\frac{1}{21} + 0.5*\frac{20}{21}より P(A|X) = 0.087となる.
これを図にすると下のようになる.
f:id:umashika5555:20170925165930p:plain
「Steveは内向的な性格だ」という事前情報を得たときの「Steveが司書である」という事後確率は少し高くなった.
コードは以下の通りである.

from IPython.core.pylabtools import figsize
import numpy as np
import matplotlib.pyplot as plt
figsize(12.5,4)
colors = ["#348ABD","#A60628"]
pr_farmer = 1/21
pr_librarian = 1 - pr_farmer
prior = [pr_farmer,pr_librarian]#事前確率
posterior = [0.087,1-0.087]#事後確率
plt.bar([0,.7],prior,alpha=0.7,width=0.25,color=colors[0],label="Prior distribution",lw="3",edgecolor="#348ABD")#事前確率
plt.bar([0+0.25,.7+0.25],posterior,alpha=0.7,width=0.25,color=colors[1],label="Posterior distribution",lw="3",edgecolor="#A60628")#事後確率
plt.xticks([0.2,0.95],["Librarian","Farmer"])#司書, 農家
plt.ylabel("Probability")#確率
plt.legend(loc="upper left")
plt.title("Prior and posterior probabilities of Steve's occupation")
plt.show()

今回は P(X|not A)を0.5に設定したがそれを0-1.0に変化させたときにどう変わっていくかを観察する.
つまり「Steveが農家だという情報が分かっているときに, Steveが内向的である確率」P(X|notA)の変化によって, Steveが司書、農家である事後確率はどのように変化するのか?

from IPython.core.pylabtools import figsize
import numpy as np
import matplotlib.pyplot as plt

figsize(12.5,4)

plt.figure()
plt.subplots_adjust(wspace=0.4, hspace=0.6)

colors = ["#348ABD","#A60628"]
pr_farmer = 1/21
pr_librarian = 1 - pr_farmer
prior = [pr_farmer,pr_librarian]#事前確率

a = np.linspace(0,1,11)
print(a)
for i,_ in enumerate(a):
    #Steveが農家であるときに内向的である事後確率
    p_x = 0.95*(1/21) + _ *(20/21)
    pos = (0.95*(1/21))/p_x
    posterior = [pos,1-pos]#事後確率
    print(i)
    plt.subplot(6,2,i+1)
    plt.bar([0,.7],prior,alpha=0.7,width=0.25,color=colors[0],label="Prior distribution",lw="3",edgecolor="#348ABD")#事前確率
    plt.bar([0+0.25,.7+0.25],posterior,alpha=0.7,width=0.25,color=colors[1],label="Posterior distribution",lw="3",edgecolor="#A60628")#事後確率

    plt.title("P(X|notA) = "+str(_))
plt.show()

f:id:umashika5555:20170925172951p:plain
これは P(A|X) = \frac{P(X|A)P(A)}{P(X)} = \frac{0.95*\frac{1}{21}}{0.95*\frac{1}{21}+x*\frac{20}{21}}=\frac{0.95}{0.95+20x}となるため,
x=0のとき「Steveが内向的だという情報を得たときのSteveが司書である」事後確率は1.0となる.
しかし,0.2以降は事後確率があまりかわらないことがわかった.
試しにグラフを書いてみた.

import numpy as np
import matplotlib.pyplot as plt
x =  np.linspace(0,1.,1000)
y = 0.95/(0.95+20*x)
plt.scatter(x,y)
plt.show()

f:id:umashika5555:20170925173414p:plain