「続・わかりやすい パターン認識」
第2章 事前確率と事後確率>2.1 事後の確率の計算>[2] コインをn回投げる場合 (p.29)

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

問題の条件はコインを1回投げる場合と同じで、 取り出したコインを複数回投げた時に表が出た回数から、取り出したのがどのコイン だったかを推定する:

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

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

  n = 10
  (0..n).each do |r|
    puts "*** #{r}"
    calc_posterior_probability_n(pi, theta, n, r).each_with_index do |post_prob, i|
      puts "P(omega_#{i+1} | H) = #{post_prob}"
    end
  end
end

# 箱の中から取り出したコインをn回投げてr回表が出た時に、
# 取り出したコインがどれだったか、それぞれのコインの確率を求める
# _pi_ :: それぞれのコインが箱の中にどのような割合で混ざっているか、の配列
# _theta_ :: それぞれのコインを投げた時に表が出る確率、の配列
# _n_ :: 試行回数
# _r_ :: 表が出た回数
def calc_posterior_probability_n(pi, theta, n, r)
  c = pi.size
  # P(x^(n)): コインをn回投げてr回表が出たという、順序ありの観測結果x^(n)が得られる確率
  head_prob = sigma(0...c) do |i|
    pi[i] * (theta[i] ** r) * ((1 - theta[i]) ** (n - r))
  end

  # それぞれのコインが取り出された確率を求める
  (0...c).map do |i|
    (pi[i] * (theta[i] ** r) * ((1 - theta[i]) ** (n - r))) / head_prob
  end
end

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

main

10回投げた結果

10回投げて表が出た回数による、取り出したコインがどれだったかの確率は、

表が出た回数 コイン1 コイン2 コイン3
0/10 0.000 0.003 0.997
1/10 0.000 0.010 0.990
2/10 0.000 0.035 0.965
3/10 0.001 0.113 0.887
4/10 0.004 0.307 0.689
5/10 0.020 0.597 0.383
6/10 0.069 0.787 0.144
7/10 0.182 0.777 0.041
8/10 0.381 0.610 0.009
9/10 0.624 0.375 0.002
10/10 0.816 0.184 0.000


graph