コインを投げるごとにどのコインなのかの確率を計算(ベイズ更新)

2015-09-02
blog

続・わかりやすい パターン認識
第2章 事前確率と事後確率>2.2 ベイズ更新 (p.32)

を、プログラムで計算してみる。

コインを複数回投げた場合では 回投げ終わって回表が出た、という具合にすべての結果が出た後で 各コインである確率を計算したが、今回は1回投げるごとに逐次更新していく(ベイズ更新)。 前回回の事後確率 を今回回目の事前確率として、 今回得られたコインの面から事後確率を更新する:
# coin_probability3.rb
class BayesUpdater
attr_reader :probs

def initialize(pi, theta)
@probs = pi.dup
@theta = theta
end

# n回目の試行がxnだった場合に、確率を更新する
def update(xn)
c = @probs.size
# n-1回目までの観測結果に対して、今回の目が出るトータルの確率
# P(x^{(n)}) = ∑_{j=1}^3 P(ω_j|x^{(n-1)}) P(x_n|ω_j)
prob = sigma(0...c) do |j|
@probs[j] * (xn * @theta[j] + (1 - xn) * (1 - @theta[j]))
end

# 各コインに対して、今回の目が出る確率を更新
c.times do |i|
@probs[i] *= (xn * @theta[i] + (1 - xn) * (1 - @theta[i])) / prob
end
end
end

def main
# πi: それぞれのコインが箱の中にどのような割合で混ざっているか
pi = [0.1, 0.4, 0.5]

# θi: それぞれのコインを投げた時に表が出る確率
theta = [0.8, 0.6, 0.3]

# 10回投げた結果が10回中7回表(HHHHTHHTHT)だった場合
#bayes_updater = BayesUpdater.new(pi, theta)
#[1, 1, 1, 1, 0, 1, 1, 0, 1, 0].each do |x|
# bayes_updater.update(x)
# p bayes_updater.probs
#end

bayes_updater = BayesUpdater.new(pi, theta)
n = 100
# 0番目のコインをn回投げてみる
t = theta[0]
head_count = 0
n.times do
x = rand < t ? 1 : 0
bayes_updater.update(x)
head_count += x
end

puts "Head: #{head_count}/#{n}"
p bayes_updater.probs
end

# Rangeに対してブロックを呼び出し、その結果の和を返す:∑
def sigma(range, &block)
range.map do |i|
block.call(i)
end.inject(:+)
end

main

# Head: 81/100
# [0.9998409538483062, 0.00015904615169388522, 3.409780749699365e-24]

10回投げて7回表が出る時の確率の遷移

コイン1 コイン2 コイン3
事前 0.1 0.4 0.5
1: 表 0.170 0.511 0.319
2: 表 0.253 0.569 0.178
3: 表 0.339 0.572 0.089
4: 表 0.423 0.535 0.042
5: 裏 0.258 0.653 0.089
6: 表 0.330 0.627 0.043
7: 表 0.404 0.576 0.020
8: 裏 0.249 0.709 0.042
9: 表 0.313 0.667 0.020
10: 裏 0.182 0.777 0.041


image