ベイジアンネットワークの例題をやり直す

以前にやったFamily Out Problemの事後確率計算は、ベイジアンネットワークらしくない計算だったことに気付いた。ネットワーク全体をまとめて計算するのでは、ネットワークの構造を無視して確率モデル全体に対してベイズの定理を使ってるだけである。ベイジアンネットワークにしてるからには、新たな情報が追加される度に、ノード毎の事前確率が事後確率に更新されるべきであろう。

問題を再掲する。


Family Out Problem

これのP(fo | lo, hb)を計算する。

この初期状態の事前確率は、単純に上から全パターンの確率を足し合わせて、
(P(fo), P(bp), P(lo), P(do), P(hb)) = (0.15, 0.01, 0.1325, 0.3958, 0.2831)
と求まる。

次に、loが確定した(lo=True, ¬lo=False)として、それに合わないP(fo)を事後確率に更新する。
ベイズの定理より、条件付き確率P(fo | lo)は
P(fo|lo) = \frac{P(lo|fo)P(fo)}{P(lo)} = \frac{P(lo|fo)P(fo)}{P(lo|fo)P(fo) + P(lo|\lnot fo)P(\lnot fo)}
= (0.15 * 0.6) / (0.15 * 0.6 + 0.85 * 0.05) = 0.6792...
であり、これに事後確率P(lo)=1.0を掛けたものがfoの事後確率なので、これを新たなP(fo)とする。
P(do), P(hb)はP(fo)に依存するので、これらも計算し直すと、
(P(fo), P(bp), P(lo), P(do), P(hb)) = (0.6792, 0.01, 1.0, 0.7103, 0.5001)
となる。

次に、hbが確定したとして、それに合わせるべく、関係する事前確率を事後確率に更新する。

まず、hbが直接依存するdoの条件付き確率は、しつこいがベイズの定理より、
P(do|hb) = \frac{P(hb|do)P(do)}{P(hb)} = \frac{P(hb|do)P(do)}{P(hb|do)P(do) + P(hb|\lnot do)P(\lnot do)}
= (0.7 * 0.7103) / (0.7 * 0.7103 + 0.01 * 0.2897) = 0.9942...
であり、hbの事後確率が1.0なので、doの事後確率がこれになる。

次に、doが依存するfoの事後確率は、¬doの事後確率 ≠ 0.0なので、¬doを考慮して P(fo|do)P(do) + P(fo|¬do)P(¬do) と計算することになる。見た目にもかなりしつこいがベイズの定理より、条件付き確率は
P(fo|do) = \frac{P(do|fo)P(fo)}{P(do)} = \frac{P(do|fo,bp)P(fo)P(bp)+P(do|fo,\lnot bp)P(fo)P(\lnot bp)}{P(do|fo,bp)P(fo)P(bp)+P(do|fo,\lnot bp)P(fo)P(\lnot bp)+P(do|\lnot fo,bp)P(\lnot fo)P(bp)+P(do|\lnot fo,\lnot bp)P(\lnot fo)P(\lnot bp)}
= (0.99 * 0.6792 * 0.01 + 0.90 * 0.6792 * 0.99) / (0.99 * 0.6792 * 0.01 + 0.90 * 0.6792 * 0.99 + 0.97 * 0.3208 * 0.01 + 0.3 * 0.3208 * 0.99) = 0.86147...
P(fo|\lnot do) = \frac{P(\lnot do|fo)P(fo)}{P(\lnot do)} = \frac{P(\lnot do|fo,bp)P(fo)P(bp)+P(\lnot do|fo,\lnot bp)P(fo)P(\lnot bp)}{P(\lnot do|fo,bp)P(fo)P(bp)+P(\lnot do|fo,\lnot bp)P(fo)P(\lnot bp)+P(\lnot do|\lnot fo,bp)P(\lnot fo)P(bp)+P(\lnot do|\lnot fo,\lnot bp)P(\lnot fo)P(\lnot bp)}
= (0.01 * 0.6792 * 0.01 + 0.10 * 0.6792 * 0.99) / (0.01 * 0.6792 * 0.01 + 0.10 * 0.6792 * 0.99 + 0.03 * 0.3208 * 0.01 + 0.7 * 0.3208 * 0.99) = 0.23232...
なので、foの事後確率はdo及び¬doの事後確率を用いて、
P(fo|do)P(do) + P(fo|¬do)P(¬do) = 0.86147 * 0.9942 + 0.23232 * (1 - 0.9942) = 0.8578...
と計算できる。

bpの事後確率も同様に計算できるので、ついでに計算してまとめると、
(P(fo), P(bp), P(lo), P(do), P(hb)) = (0.8579, 0.0138, 1.0, 0.9942, 1.0)
となる。

検算のために書いたPythonのコード 

#prior prob
P1 = 0.15
P2 = 0.01

def P_fo(fo):
    return (1 - P1, P1)[fo]

def P_bp(bp):
    return (1 - P2, P2)[bp]

def P_lo(lo, fo):
    if fo != None:
        return ((0.95, 0.05), (0.4, 0.6))[fo][lo]
    else:
        return sum([P_lo(lo, fo) * P_fo(fo) for fo in [True, False]])

def P_do(do, fo, bp):
    if fo != None and bp != None:
        return (((0.7, 0.3), (0.03, 0.97)), ((0.1, 0.9), (0.01, 0.99)))[fo][bp][do]
    elif fo != None and bp == None:
        return sum([P_do(do, fo, bp) * P_bp(bp) for bp in [True, False]])
    elif fo == None and bp != None:
        return sum([P_do(do, fo, bp) * P_fo(fo) for fo in [True, False]])
    else:
        return sum([P_do(do, fo, bp) * P_fo(fo) * P_bp(bp) for fo in [True, False] for bp in [True, False]])

def P_hb(hb, do):
    if do != None:
        return ((0.99, 0.01), (0.3, 0.7))[do][hb]
    else:
        return sum([P_hb(hb, do) * P_do(do, None, None) for do in [True, False]])

#prior prob
P3 = P_lo(True, None)
P4 = P_do(True, None, None)
P5 = P_hb(True, None)
print("P(fo)={:.4f}, P(bp)={:.4f}, P(lo)={:.4f}, P(do)={:.4f}, P(hb)={:.4f}".format(P1, P2, P3, P4, P5))

#posterior prob given P3=1.0
P3 = 1.0
P1 = P_lo(True, fo=True) * P_fo(True) / P_lo(True, None)
P4 = P_do(True, None, None)
P5 = P_hb(True, None)
print("P(fo)={:.4f}, P(bp)={:.4f}, P(lo)={:.4f}, P(do)={:.4f}, P(hb)={:.4f}".format(P1, P2, P3, P4, P5))

#posterior prob given P5=1.0
P5 = 1.0
P4 = P_hb(True, do=True) * P_do(True, None, None) / P_hb(True, None) 
#P1 = P(fo=True|do=True)P(do=True) + P(fo=True|do=False)P(do=False)
P1, P2 = P_do(True, fo=True, bp=None) * P_fo(True) / P_do(True, None, None) * P4 \
       + P_do(False, fo=True, bp=None) * P_fo(True) / P_do(False, None, None) * (1 - P4), \
         P_do(True, fo=None, bp=True) * P_bp(True) / P_do(True, None, None) * P4 \
       + P_do(False, fo=None, bp=True) * P_bp(True) / P_do(False, None, None) * (1 - P4)
print("P(fo)={:.4f}, P(bp)={:.4f}, P(lo)={:.4f}, P(do)={:.4f}, P(hb)={:.4f}".format(P1, P2, P3, P4, P5))

出力 
P(fo)=0.1500, P(bp)=0.0100, P(lo)=0.1325, P(do)=0.3958, P(hb)=0.2831
P(fo)=0.6792, P(bp)=0.0100, P(lo)=1.0000, P(do)=0.7103, P(hb)=0.5001
P(fo)=0.8579, P(bp)=0.0138, P(lo)=1.0000, P(do)=0.9942, P(hb)=1.0000