Estimación de la tasa de fallos a partir de los datos observados

Question

Hace poco leí el excelente libro Programación probabilística y métodos bayesianos para hackers y estoy tratando de resolver algunos problemas por mi cuenta:

Realizo un experimento para estimar la fiabilidad de un dispositivo (por ejemplo CPU, coche, ventilador, etc.) comprando 1000 dispositivos y haciéndolos funcionar durante un período de 1 año para ver cuántos fallan. Mi experimento muestra que 0 de de ellos han fallado durante ese periodo. ¿Cuál es mi intervalo de confianza del 95% de la tasa de fallos?

Si el número de dispositivos que fallaron fuera distinto de cero, podría estimar la media y el intervalo de confianza con técnicas estándar como el cálculo del error estándar o el bootstrap. El hecho de que hayan fallado 0 dispositivos lo hace mucho más complicado. Intento utilizar PyMC para resolverlo:

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. ¿Es correcto mi modelo?
2. ¿Es correcto mi uso de PyMC?
3. ¿Puede alguien confirmar mi respuesta con otro método? ¿Tal vez de forma analítica utilizando una distribución de Poisson?

P.D. Prácticamente no sé nada de la teoría que hay detrás de MCMC, pero el enfoque de la aplicación, primero, y de las matemáticas, después, del libro es excelente.

Answer 1

Este parece un buen modelo, implementado correctamente en PyMC. Hay dos hechos estadísticos bayesianos que podemos utilizar para confirmar su respuesta con otro método:

1. $\textrm{Beta}(1,1)$ es equivalente a la distribución uniforme en el intervalo $[0,1]$ ;
2. Las distribuciones beta y binomial son conjugadas.

Esto significa que la distribución posterior de $p$ es también una distribución beta, y (si he entendido bien la parametrización) $p_{\text{posterior}} \sim \textrm{Beta}(1,1001)$ . Puede comparar los percentiles de esta distribución derivada analíticamente con los percentiles que ha encontrado mediante MCMC de esta manera:

> 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]

Aquí hay un cuaderno que recoge todo esto .

Por cierto, no estoy seguro de que los estadísticos con carné llamen a esto un intervalo de confianza. Hay muchas sutilezas en este asunto, pero la principal objeción práctica es que no has especificado cómo los 1000 dispositivos fueron seleccionados. Por ejemplo, si usted es un comprador muy importante, y está probando estos mil antes de pedir un millón de unidades, el fabricante podría hacer todo lo posible para asegurarse de que obtiene un buen millar en este lote.