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()