bayes.py

import scipy.stats as stats
import numpy as np
import matplotlib.pyplot as plt

# create model and obtain posterior
n_heads = 6
n = 9
grid_size = 1000
p = np.linspace(0, 1, grid_size)
prior = stats.beta.pdf(p, 10, 10)
likelihood = stats.binom.pmf(n_heads, n, p)
posterior = likelihood * prior
posterior = posterior / np.sum(posterior)


ff, axarr = plt.subplots(3, sharex=False)
ff.tight_layout()
ff.subplots_adjust(left=0.1)
axarr[0].plot(p, prior)
axarr[0].set_title('Prior')
axarr[1].plot(p, likelihood)
axarr[1].set_title('Likelihood')
axarr[2].plot(p, posterior)
axarr[2].set_title('Posterior')
ff.savefig('dists.png')

# sampling from posterior distribution
sample_size = 1e4
samples = np.random.choice(p, size=sample_size, replace=True, p=posterior)

# Calculate quantities of interest
mu = samples.mean()
var = samples.var()
lower, upper = np.percentile(samples, [5, 95])

print mu
print var
print lower, upper

# estimate the density
density = stats.kde.gaussian_kde(samples)

ff, axarr = plt.subplots(2, sharex=False)
axarr[0].scatter(np.arange(sample_size), samples)
axarr[0].set_title('Samples from posterior')
axarr[0].set_xlim(0, sample_size)
axarr[1].plot(p, density(p))
axarr[1].set_title('Density of samples from posterior')
axarr[1].fill_between(p, density(p), color='#3399ff', where=(p > lower) & (p < upper))
ff.savefig('samples.png')