In [1]:
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)

In [2]:
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
In [3]:
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]


# Analytic Solution: beta dist and binomial dist are conjugate¶

In [4]:
# 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)

In [5]:
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)')

Out[5]:
<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]))

In [6]:
# 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)

In [7]:
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)')

Out[7]:
<matplotlib.axes.AxesSubplot at 0x60dea10>
In [8]:
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]

In [8]: