写経

import random
import pylab

pylab.rc('figure', dpi=200)

20.3 ベイズ更新

code 20.1

def calc_bays(prior_A, prob_b_if_A, prob_B):
    return prior_A*prob_b_if_A/prob_B

prior_A = 1/3
prob_6_if_A = 1/5
prob_6 = (1/5 + 1/6 + 1/7)/3

post_A = calc_bays(prior_A, prob_6_if_A, prob_6)
print(f'タイプAである確率 = {round(post_A, 4)}')
post_A = calc_bays(post_A, prob_6_if_A, prob_6)
print(f'タイプAである確率 = {round(post_A, 4)}')

タイプAである確率 = 0.3925
タイプAである確率 = 0.4622

post_A = calc_bays(prior_A, 1 - prob_6_if_A, 1 - prob_6)
print(f'タイプAである確率 = {round(post_A, 4)}')
post_A = calc_bays(post_A, 1 - prob_6_if_A, 1 - prob_6)
print(f'タイプAである確率 = {round(post_A, 4)}')

タイプAである確率 = 0.3212
タイプAである確率 = 0.3096

prior_A = 9/10
post_A = calc_bays(prior_A, 1 - prob_6_if_A, 1 - prob_6)
print(f'タイプAである確率 = {round(post_A, 4)}')
post_A = calc_bays(post_A, 1 - prob_6_if_A, 1 - prob_6)
print(f'タイプAである確率 = {round(post_A, 4)}')

タイプAである確率 = 0.8673
タイプAである確率 = 0.8358

code 20.2

def bays_update(num_rolls=200, prior_A=9/10):
    post_A = prior_A
    for i in range(num_rolls + 1):
        if i%(num_rolls//10) == 0:
            print(f'{i}回振った時点でタイプAである確率 = {round(post_A, 4)}')
        is_six = random.random() <= 1/7
        if is_six:
            post_A = calc_bays(post_A, prob_6_if_A, prob_6)
        else:
            post_A = calc_bays(post_A, 1 - prob_6_if_A, 1 - prob_6)

bays_update()

0回振った時点でタイプAである確率 = 0.9
20回振った時点でタイプAである確率 = 0.7834
40回振った時点でタイプAである確率 = 0.4567
60回振った時点でタイプAである確率 = 0.5936
80回振った時点でタイプAである確率 = 0.6314
100回振った時点でタイプAである確率 = 1.0028
120回振った時点でタイプAである確率 = 1.0667
140回振った時点でタイプAである確率 = 0.6218
160回振った時点でタイプAである確率 = 0.6614
180回振った時点でタイプAである確率 = 0.7035
200回振った時点でタイプAである確率 = 0.7483

メモ

　どうやら、事後確率が1以上になることがあるようだ。というか、そもそも

print(f'prob_6_if_A/prob_6 = {prob_6_if_A/prob_6}')

prob_6_if_A/prob_6 = 1.177570093457944

なので、random.random() <= 1/7でTrueが続くと事後確率が1以上になるのは当たり前だ。確立を論じている章としていかがなものか。

　次の章読んでこういうコードが続くならこの本は切ろうかな・・・

メモ

　ベイズ更新の様子を可視化するためにプロットを行ってみた。上でも述べたが、このベイズ更新の式は確率論的に誤りがあるので、参考程度にプロットしたまでである。

def bays_update(num_rolls=200, prior_A=9/10):
    prob_6_if_A = 1/5
    prob_6 = (1/5 + 1/6 + 1/7)/3
    post_A = prior_A
    post_As = []
    post_As.append(post_A)
    for i in range(num_rolls + 1):
        is_six = random.random() <= 1/7
        if is_six:
            post_A = calc_bays(post_A, prob_6_if_A, prob_6)
            post_As.append(post_A)
        else:
            post_A = calc_bays(post_A, 1 - prob_6_if_A, 1 - prob_6)
            post_As.append(post_A)
    return post_As

def plot_bays_update(num_rolls=200, prior_A=9/10, num_trials=10):
    for i in range(1, 1 + num_trials):
        post_As = bays_update(num_rolls, prior_A)
        pylab.plot(post_As, label=f'trial {i}')
    pylab.title('Bays uppdate')
    pylab.xlabel('Dice throw')
    pylab.ylabel('Post probability')
    pylab.legend(loc='best')