python實作 beta-binomial 的共轭分布

下圖可見, 随着N增加, 分布的均值趨于準确, 每張圖的藍色,紅色和綠色分别代表α和β為(1, 1), (0.5, 0.5), (20, 20)的先驗機率, 可見到最右下的一張圖, 先驗機率依然有可見的而影響.

# -*- coding: utf-8 -*-
import matplotlib.pyplot as plt
import numpy as np
from scipy import stats
from matplotlib import style
style.use('ggplot')

# p( theta | y ) = Beta( a + y, b + N - y )
theta_real = 0.35

#為公式中的N
trials = [0, 1, 2, 3, 4, 8, 16, 32, 50, 150]

#為公式中的y
data = [0, 1, 1, 1, 1, 4, 6, 9, 13, 48]

#每個tuple裡分别為alpha和beta的先驗機率
beta_params = [(1, 1), (0.5, 0.5), (20, 20)]

x = np.linspace(0, 1, 100)
# 'b', 'r', 'g' 分别對應 (1, 1), (0.5, 0.5), (20, 20)
colors = ['b','r','g']
for i, N in enumerate(trials):
    if i == 0:
        plt.subplot(432)
    else:
        plt.subplot(4, 3, i+3)
    y = data[i]
    c = 0
    for a_prior, b_prior in beta_params:
        p_theta_given_y = stats.beta.pdf(x, a_prior+y, b_prior+N-y)
        plt.plot(x, p_theta_given_y, colors[c])
        plt.fill_between(x, 0, p_theta_given_y, color=colors[c], alpha=0.6)
        c += 1
    plt.axvline(theta_real, ymax=0.3, color='k')
    plt.plot(0,0,label='{:d} experiments\n{:d} heads'.format(N,y),alpha=0)
    plt.xlim(0,1)
    plt.ylim(0,12)
    plt.xlabel(r'$\theta$')
    plt.legend(fontsize=8)
plt.suptitle('Beta prior with Binomial likelihood', fontsize=16)
plt.show()