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]



*

*Is my model correct?

*Is my usage of PyMC correct?

*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.
 A: This does seem like a good model, implemented correctly in PyMC.  There are two Bayesian stats facts that we can use to confirm your answer with another method:


*

*$\textrm{Beta}(1,1)$ is equivalent to the uniform distribution on the interval $[0,1]$;

*The beta and binomial distributions are conjugate.


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:
> p_posterior = np.random.beta(a=1, b=1001, size=1000000)
> print np.percentile(p_posterior, [2.5, 97.5])
[2.5350975458273468e-05, 0.003681783314872197]

Here is a notebook that collects this all up.
By the way, I'm not sure if card-carrying statisticians would call this a confidence interval.  There is a lot of subtlety to this issue, but the major practical objection is that you didn't specify how the 1000 devices were selected.  For example, if you are a very important purchaser, and you are testing these thousand before ordering a million units, the manufacturer could go to great lengths to ensure you get a good thousand in this batch!  
