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