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

# 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]