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]

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!


Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

Not the answer you're looking for? Browse other questions tagged or ask your own question.