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)
[-----------------100%-----------------] 40000 of 40000 complete in 3.1 sec
print np.percentile(mcmc.trace('p')[:], [2.5, 97.5])
print p.stats()['95% HPD interval'] # alternative approach to interval prediction
[2.6165440508793241e-05, 0.0036806413183323641] [ 7.48379837e-08 2.95055578e-03]
# 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)')
<matplotlib.axes.AxesSubplot at 0x5a9e950>
/homes/abie/anaconda/lib/python2.7/site-packages/matplotlib/font_manager.py:1236: UserWarning: findfont: Font family ['Arial'] not found. Falling back to Bitstream Vera Sans (prop.get_family(), self.defaultFamily[fontext]))
# 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)')
<matplotlib.axes.AxesSubplot at 0x60dea10>
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
From MCMC [ 3.00000000e-05 3.68000000e-03] [ 0. 0.00295] From Beta(1,1001) [ 3.00000000e-05 3.68000000e-03] [ 0. 0.00299]