# 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.

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:

1. $\textrm{Beta}(1,1)$ is equivalent to the uniform distribution on the interval $[0,1]$;
2. 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]


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!