¿Cómo estimar la distribución de Poisson con una ocurrencia de eventos?

Question

Digamos que observamos un evento el 1 de enero de 2005, y que observamos los eventos desde el año 2000. ¿Cuál es la mejor estimación de la distribución de Poisson?

Un enfoque es establecer la intensidad $\lambda=1/18$ al año. Se trata de una estimación MLE de Poisson, ya que han pasado 18 años y sólo se ha observado un evento.

Otro enfoque parece ser el ajuste de la distribución exponencial al tiempo de llegada. ¿Conducen a la misma estimación?

ACTUALIZACIÓN

En mi problema la estimación de la intensidad de Poisson está condicionada a la observación de un único evento. Este es un factor importante, porque no es obvio que no estropee el MLE habitual para Poisson.

He realizado una simulación en la que se extraen intensidades aleatorias, y luego se extraen frecuencias aleatorias a partir de estas intensidades. Sólo miramos las frecuencias iguales a uno de cada 18 años. Finalmente, utilizamos el estimador de intensidad que es igual a la frecuencia observada, como en un MLE habitual para Poisson, y lo comparamos con la intensidad verdadera.

Sospeché que para los eventos raros, cuando el evento se observó precisamente una vez en un período de la muestra, podría haber un mejor estimador de la intensidad de Poisson que la frecuencia observada. Sin embargo, esta simulación muestra un sesgo nulo de la MLE de Poisson habitual.

import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as stats

ns = 100000
T = 10
lam = np.random.gamma(1,1/T,ns) # random intensities ~1/18 (per year)
obsf = np.random.poisson(lam*T,ns) # observe frequencies for 18 years

plt.hist(obsf)
plt.title('random intensities')
plt.show()

sidx = np.where(obsf==1)[0] # we only look at trials where 1 event happened in 18 years
err = (obsf[sidx]/T - lam[sidx]) # error when intensity estimator is set to observed frequency
print(stats.describe(err))

plt.hist(err)
plt.title('errors')
plt.show()

de salida:

DescribeResult(nobs=25146, minmax=(-0.30917753984994045, 0.055486021558327742), mean=-0.00044121203771743959, variance=0.0015661042110224508, skewness=-1.434031238643681, kurtosis=2.993433694497429)

ACTUALIZACIÓN 2

Utilicé Gamma como la distribución de parámetros simulada. Es injusto, porque Gama es una previa de Poisson. Usaré alguna otra distribución.

Answer 1

Yo trataría esto como una observación de, suponiendo que se empezó a medianoche del 31 de diciembre de 2000, 4 años, y otra observación censurada a 12,956 años. En el caso de la distribución exponencial, este es un problema bastante fácil de resolver:

$$p(x_1, x_2 | \lambda) = \left(\frac{1}{\lambda}e^{-x_1/\lambda}\right)\left(e^{-c_2/\lambda}\right)$$

donde el primer término es la probabilidad de observar $x_1$ y la segunda la probabilidad de observar $x_2 > c_2$ con $c_2$ siendo el punto de censura de 12,956 años. Tomando el logaritmo y poniendo la derivada igual a cero nos da:

$$\hat{\lambda} = \frac{1}{x_1 + c_2}$$

lo que nos da aproximadamente $1/18$ años, ya que estamos a finales de 2017 en vez de a principios.

Esto es sólo un caso específico del resultado general de que la MLE para una distribución exponencial censurada es igual a la suma de los tiempos totales observados, independientemente de si fueron censurados, dividida por el número total de llegadas observadas.

La conclusión: los dos enfoques darán el mismo resultado.

Editar en respuesta a los comentarios:

Los estimadores basados en Poisson y Exponencial siempre darán la misma respuesta. En el caso del Poisson, se tiene $n$ eventos observados durante el periodo de tiempo $T$ dando una estimación de $T/n$ . En el caso del Exponencial, su tiempo total observado es $T$ y tienes $n+1$ observaciones, con la última censurada al final del periodo de observación, por lo que la estimación es $n/T$ .

Answer 2

Mirando más de cerca tu simulación, puedo ver un poco más claramente lo que estás haciendo. Has configurado un proceso de generación de datos con un intensidad aleatoria tomando una distribución gamma para cada intensidad gamma, hay una realización de Poisson, entonces se pasa a estimar la tasa de Poisson entre una submuestra con recuentos no nulos.

Esto es un poco diferente de cómo me imaginaba el problema de las fechas porque tienes datos de una muestra, algunos de los cuales están deliberadamente excluidos, y esta muestra tiene una importante heterogeneidad.

Hay un par de cosas que abordar aquí: Si no se tiene una entrada aleatoria al proceso de Poisson, el modelo de probabilidad correcto para los datos observados es un distribución de Poisson truncada por cero . https://en.wikipedia.org/wiki/Zero-truncated_Poisson_distribution . Obsérvese que un estimador del método de los momentos para la intensidad podría resolverse numéricamente (no hay solución analítica manejable)

\begin{equation} 0 = \bar{Y} - \lambda \exp(\lambda) / (\exp(\lambda) - 1) \end{equation}

set.seed(123)
ns <- 10000
zcpois <- function(y) uniroot(function(x) mean(y) - x*exp(x)/(exp(x)-1), c(1e-8, max(y)))
y <- rpois(ns, 0.7)
zcpois(y[y!=0]) ## very precise MoM

> zcpois(y[y!=0])$root
[1] 0.6950045

También se puede utilizar la máxima verosimilitud recordando la verosimilitud para un VR truncado: $P(Y=y | Y>0) = P(Y=y)/P(Y>0)$ .

y0 <- y[y!=0]
negloglik <- function(lambda, y) {
  -sum(dpois(x=y, lambda=lambda, log=T) - ppois(q=0, lambda=lambda, lower.tail=F, log=T))
}

lseq <- seq(from=0.001, to=2, by=0.001)
nll <- sapply(lseq, negloglik, y=y0)
plot(lseq, nll, type='l')
i <- which.min(nll)
abline(v=lseq[i])
text(lseq[i], max(nll), paste('MLE:', lseq[i]), pos=4, xpd=T)

Sin embargo, es importante considerar la mezcla Poisson/Gamma y el impacto en el modelo de probabilidad resultante. El modelo de mezcla no sigue una distribución de Poisson truncada en cero. El proceso de Poisson subyacente (sin omitir los recuentos 0) seguiría un proceso de Poisson sobredispersado.

lambdas <- rgamma(ns, 1, 18) ## rate parametrization
counts <- rpois(ns, lambdas)
nzcounts <- counts[counts!=0]

> zcpois(nzcounts)$root
[1] 0.1311875
> mean(lambdas)
[1] 0.05573894

Un modelo bayesiano cuya verosimilitud sigue un modelo de Poisson truncado en cero y cuya intensidad sigue una distribución gamma con un parámetro de forma y escala no informativos le permitiría estimar la distribución posterior, cuyos modos posteriores le darían buenas estimaciones de los parámetros de entrada a la distribución gamma. Un enfoque frecuentista invocaría el algoritmo EM.