ベイズ誤り確率の計算

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

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

  n = 10
  eb = sigma(0..n) do |r|
    # コインをn回投げて表がr回出た時に、投げたコインが1~3であった確率
    probs = calc_posterior_probability_n(pi, theta, n, r)

    # ベイズ決定則で、確率が最大のものを選択する
    max = probs.max

    # コイン1~3をn回投げた時に表がr回出る確率
    pnr = sigma(0...pi.size) do |i|
      bernoulli_distribution(n, r, theta[i]) * pi[i]
    end

    #puts "*** #{r}"
    #puts "#{r}: #{max}, #{1 - max}, #{pnr}"

    # r回表が出る時に判定を間違う確率と、r回表が出るということが起こる確率を掛ける
    ebr = (1 - max) * pnr
  end
  puts eb
end

# 2.2 より
def calc_posterior_probability_n(pi, theta, n, r)
  m = pi.size
  head_prob = (0...m).map do |i|
    pi[i] * (theta[i] ** r) * ((1 - theta[i]) ** (n - r))
  end.inject(:+)

  (0...m).map do |i|
    (pi[i] * (theta[i] ** r) * ((1 - theta[i]) ** (n - r))) / head_prob
  end
end

# ベルヌーイ分布
def bernoulli_distribution(n, r, theta)
  combi(n, r) * (theta ** r) * ((1 - theta) ** (n - r))
end

# 組み合わせ nCr
def combi(n, r)
  fact(n) / (fact(n - r) * fact(r))
end

# 階乗 n!
def fact(n)
  return 1 if n <= 1
  (1..n).inject(:*)
end

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

main

# n=10   => 0.22259161014
# n=100  => 0.006150932942506893