我最近读了一本优秀的书Probabilistic Programming & Bayesian Methods for Hackers，我正在尝试自己解决一些问题：
我通过购买 1000 台设备并在 1 年内运行它们以查看有多少失败，来进行一个实验来估计设备（例如 CPU、汽车、风扇等）的可靠性。我的实验表明，其中有 0 个在此期间失败了。失败率的 95% 置信区间是多少？
如果失败的设备数量不为零，我可以使用计算标准误差或自举等标准技术来估计平均值和置信区间。0 个设备失败的事实使它变得更加棘手。我尝试使用 PyMC 来解决它：
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]
- 我的模型正确吗？
- 我对 PyMC 的使用正确吗？
- 有人可以用另一种方法确认我的答案吗？也许通过使用泊松分布进行分析？
PS 我实际上对 MCMC 背后的理论一无所知，但是这本书的应用优先和数学第二的方法非常好。