続・わかりやすい パターン認識
第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