2021/01/11
今回は
前回、
前々回に手計算した内容をRで実装します。Pythonでも良かったのですが、Rの方が参考資料が多く、手っ取り早く実装できそうだったので選びました。
実装したRのコードは以下の通りです。
library(bnlearn)
# Bayesian network
dag = model2network("[C][S|C][R|C][W|S:R]")
plot(dag)
# Probability
C.lv = c("F", "T")
C.prob = array(c(0.50, 0.50), dim = 2, dimnames = list(C = C.lv))
S.lv = c("F", "T")
S.prob = array(c(0.50, 0.50, 0.90, 0.10), dim = c(2, 2), dimnames = list(S = S.lv, C = C.lv))
R.lv = c("F", "T")
R.prob = array(c(0.80, 0.20, 0.20, 0.80), dim = c(2, 2), dimnames = list(R = R.lv, C = C.lv))
W.lv = c("F", "T")
W.prob = array(c(1.00, 0.00, 0.10, 0.90, 0.10, 0.90, 0.01, 0.99), dim = c(2, 2, 2), dimnames = list(W = W.lv, S = S.lv, R = R.lv))
cpt = list(C = C.prob, S = S.prob, R = R.prob, W = W.prob)
# Model fit
survey = custom.fit(dag, cpt)
# Prediction
N = 2000
temp = array(0,N)
# P(C=T|W=T)
for (i in 1:N) {
temp[i] = cpquery(survey, event = (C == "T"), evidence = (W == "T"))
}
print(paste("P(C=T|W=T) =", mean(temp)))
# P(C=T|R=T)
for (i in 1:N) {
temp[i] = cpquery(survey, event = (C == "T"), evidence = (R == "T"))
}
print(paste("P(C=T|R=T) =", mean(temp)))
# P(S=T|R=T)
for (i in 1:N) {
temp[i] = cpquery(survey, event = (S == "T"), evidence = (R == "T"))
}
print(paste("P(S=T|R=T) =", mean(temp)))
# P(W=T|R=T)
for (i in 1:N) {
temp[i] = cpquery(survey, event = (W == "T"), evidence = (R == "T"))
}
print(paste("P(W=T|R=T) =", mean(temp)))
実行結果は以下の通り。
5行目のプロットはこんな感じになりました。線がクロスしていますが、前回、前々回と議論してきたモデルです。

8〜16行目は、確率の値を代入し、19行目でモデルフィットを行っています。22行目以降は手計算で行った条件付き確率を求めています。計算はcpquery関数を使っていますが、この関数は以下の注意点がマニュアルにありましたので、2,000回計算させて、平均値を求めました。
「cpqueryとcpdistの両方ともモンテカルロ粒子フィルタに基づいているため,シミュレーションのノイズの影響で,実行によってはわずかに異なる値を返す可能性があることに注意してください.」
P(C=T|W=T) = 0.575912086190798(手計算:0.5758)
P(C=T|R=T) = 0.799937096528357(手計算:0.8)
P(S=T|R=T) = 0.179835101789838(手計算:0.18)
P(W=T|R=T) = 0.916275204803596(手計算:0.9162)
簡単に条件付き確率を求めることができたので、拍子抜けした感がありますが、変数(ノード)が増えると、さすがに手計算では無理ですね。Rで簡単に実装できることが分かりましたので、何か別のテーマを検討してみようと思います。