# Estimating Failure Rate from Observed Data

I recently read the excellent book Probabilistic Programming & Bayesian Methods for Hackers and I'm trying to solve some problems on my own:

I perform an experiment to estimate the reliability of a device (e.g. CPU, car, fan, etc.) by purchasing 1000 devices and running them over a 1 yr period to see how many fail. My experiment shows that 0 of them have failed over the period. What is my 95% confidence interval of the failure rate?

If the number of devices that failed were non-zero, I could estimate the mean and confidence interval with standard techniques like calculating the standard error or bootstrapping. The fact that 0 devices failed makes it much trickier. I try and use PyMC to solve it:

import numpy as np
import pymc as pm

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

> [3.3206054853225512e-05, 0.0037895137242935613]
1. Is my model correct?
2. Is my usage of PyMC correct?
3. Can someone confirm my answer with another method? Maybe analytically by using a Poisson distribution?

P.S. I practically know nothing about the theory behind MCMC but the application-first and mathematics-second approach of the book is excellent.

1. $\textrm{Beta}(1,1)$ is equivalent to the uniform distribution on the interval $[0,1]$;
This means that the posterior distribution of $p$ is also a beta distribution, and (if I have got the parameterization correct) $p_{\text{posterior}} \sim \textrm{Beta}(1,1001)$. You can compare the percentiles from this analytically derived distribution with the percentiles that you have found via MCMC thusly: