import numpy as np, pymc as pm, seaborn as sns
%matplotlib inline
import matplotlib as mpl, matplotlib.pyplot as plt, mpld3
mpld3.enable_notebook()

mpl.rcParams['figure.figsize'] = (10, 4)
np.random.seed(123456)

data = np.zeros(1000)  # observed data: zero failures out of 1000 devices
p = pm.Uniform('p', 0, 1) # model the failure rate as a uniform distribution from 0 to 1
obs = pm.Bernoulli('obs', p, value=data, observed=True) # each device can fail or not fail. i.e. the observations follow a Bernoulli distribution
model = pm.Model([obs, p])

mcmc = pm.MCMC(model)
mcmc.sample(40000, 10000, 1)

print np.percentile(mcmc.trace('p')[:], [2.5, 97.5])
print p.stats()['95% HPD interval']  # alternative approach to interval prediction

# Uniform(0,1) is same distribution at Beta(1,1)

u = np.random.uniform(low=0, high=1, size=1000000)
b = np.random.beta(a=1, b=1, size=1000000)

fig, ax = plt.subplots(1, 2, sharex=True, sharey=True)
sns.distplot(u, bins=10, ax=ax[0], kde=True, axlabel='Uniform(0,1)')
sns.distplot(b, bins=10, ax=ax[1], kde=True, axlabel='Beta(1,1)')

# After observing 0 successes from 1000 binomially distributed trials, the posterior of p is still beta-distributed
# with p_posterior ~ Beta(1,1001)

p_posterior = np.random.beta(a=1, b=1001, size=1000000)

fig, ax = plt.subplots(1, 2, sharex=True, sharey=True)
sns.distplot(p.trace(), bins=10, ax=ax[0], kde=True, axlabel='MCMC Samples')
sns.distplot(p_posterior, bins=10, ax=ax[1], kde=True, axlabel='Beta(1,1001)')

print 'From MCMC'
print np.round(np.percentile(mcmc.trace('p')[:], [2.5, 97.5]), 5)
print np.round(p.stats()['95% HPD interval'], 5)  # alternative approach to interval prediction
print
print 'From Beta(1,1001)'
print np.round(np.percentile(p_posterior, [2.5, 97.5]), 5)
print np.round(pm.utils.hpd(p_posterior, alpha=.05), 5)  # alternative approach to interval prediction