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.