La simulación de un Poisson mixto me da menos eventos de los esperados.

Question

Haciendo referencia a la sección 5.4.3 del libro Introducción a los modelos de probabilidad de Sheldon Ross (10ª edición), se introduce la Poisson mixta en la que la tasa, $\lambda$ de un proceso de Poisson se extrae a su vez de una distribución Gamma. Dejemos que esta distribución Gamma que estamos dibujando $\lambda$ de ser descrito por $L$ . En la página 352 tenemos:

$$E(N(t)|L) = Lt$$

Y así obtenemos el resultado intuitivo:

$$E(N(t)) = E(E(N(t)|L)) = E(L)t$$

Estoy tratando de replicar esto con la simulación. Empecemos con un simple proceso de Poisson. Usamos Python para simular distribuciones exponenciales entre llegadas y nos detenemos cuando ha transcurrido un tiempo total de 1000 unidades. Cuando la tasa, $\lambda=5$ Según el código siguiente, obtenemos unos 5000 eventos una vez que han transcurrido 1000 unidades de tiempo. No hay sorpresas aquí.

import numpy as np
period=1000
evnts = 0
t=0
lm = 5
while t<period:
    t_del = np.random.exponential(1/lm)
    t+=t_del
    evnts+=1
print(evnts)

Ahora, quiero extender esta simulación a una Poisson mixta como se describe en Ross. En lugar de simular a partir de una exponencial fija cada vez, primero simulo la $\lambda$ de una distribución Gamma. He confirmado que la Gamma utilizada tiene una media de 5. A continuación, utilizo esa $\lambda$ para simular una distribución exponencial. Y el número de eventos en 1000 unidades de tiempo baja a unos 3000 en lugar de 5000 como se esperaba de la ecuación anterior. ¿Qué es lo que me falta? ¿Por qué el número de eventos debería ser menor?

import numpy as np
period=1000
evnts = 0
t=0
while t<period:
    theta = 0.5
    lm = np.random.gamma(5*theta,1/theta)
    t_del = np.random.exponential(1/lm)
    t+=t_del
    evnts+=1
print(evnts)

Answer 1

Para simular un Poisson mixto se genera primero un solo Variable gamma $L$ . Este valor para $L$ le da la $\lambda$ (tasa), que se utiliza para la duración de su proceso de Poisson. En otras palabras, tienes que mover la línea

lm = np.random.gamma(5*theta,1/theta)

fuera de su while bucle. Al hacer esto, el recuento de eventos acabará acercándose a los 5000 previstos.

Su código actual no está simulando un proceso de Poisson, ya que los tiempos entre llegadas no se extraen de una única distribución exponencial. Por el contrario, cada tiempo entre llegadas se extrae de una diferentes exponencial, cuyo parámetro de tasa $\lambda$ oscila entre $5-\sqrt{10}$ a $5+\sqrt{10}$ (estos valores son la media $\pm$ desviación estándar). Su código, tal como está escrito, generará muchos valores bajos para $\lambda$ (que corresponden a tiempos de interllegada elevados, que consumen gran parte del periodo) en prácticamente todas las ejecuciones de su proceso puntual simulado; esto es lo que hace que el recuento de eventos disminuya. Cuando se simula correctamente, la baja $\lambda$ proceso se producirá, pero con mucha menos frecuencia, por lo que los recuentos de eventos bajos son mucho menos frecuentes.