Simulación de Monte Carlo usando distribuciones exponenciales

Question

Estoy tratando de simular un modelo estocástico de crecimiento exponencial determinista de la población, donde $dN/dt = rN$ donde $N$ es el tamaño de la población y $r$ es la tasa ($t$ tiempo). Estoy asumiendo que no hay capacidad de carga. Esta página (http://cnr.lwlss.net/DiscreteStochasticLogistic/) sugiere este algoritmo para simular el crecimiento en un intervalo $[0, t_{end}]$:

1. empieza en $t = 0$ con tamaño de población inicial
2. elige el siguiente tiempo para el evento de nacimiento, $\delta t \sim Exponencial(rN(1 - N/K))$ ($K$ es la capacidad de carga)
3. incrementa el tamaño de la población, $N = N + 1$
4. establece $t = t + \delta t$
5. si $t > t_{end}$ entonces termina, de lo contrario ve al paso 2.

Dado que no tengo una capacidad de carga $K$, asumo que es infinita, entonces el próximo tiempo de nacimiento es $\delta t \sim Exponencial(rN)$. ¿Es eso correcto?

cuando ejecuto la simulación de esta manera, no proporciona en absoluto resultados similares a $N(t) = P_0 e^{rt}$ (donde $P_0$ es el tamaño inicial de la población). incluso promediando sobre muchas iteraciones parece dar curvas de crecimiento diferentes.

a continuación se muestra mi código y el resultado de la simulación. la curva roja es crecimiento exponencial determinista, las curvas negras son simulaciones usando distribución exponencial. claramente no coinciden.

import numpy as np
import matplotlib.pylab as plt
def sim(rate, start, end, init):
    N = 200
    tamañosfinales = []
    resultados = []
    for n in range(N):
        tamaño = init
        curr_t = 0
        tiempos = [curr_t]
        tamaños = [init]
        nueva_tasa = rate
        while curr_t <= end:
            # simular próximo tiempo de nacimiento. el parámetro de escala
            # es inversamente proporcional al tamaño de la población
            nueva_tasa = 1/float(nueva_tasa * tamaño)
            tiempo_div = np.random.exponential(scale=nueva_tasa)
            # avanzar en el tiempo
            curr_t += tiempo_div
            if curr_t > end:
                # si excedemos el intervalo de tiempo, terminamos
                break
            tiempos.append(curr_t)
            # aumentar el tamaño de la población
            tamaño += 1
            tamaños.append(tamaño)
        tamañosfinales.append([tiempos, tamaños])
    return tamañosfinales

# ejecutar simulación y graficar resultados
init = 20
start = 0
end = 20
rate = 1
tamañosfinales = sim(rate, start, end, init)
plt.figure()
todoslostamaños = []
for f in tamañosfinales:
    todoslostamaños.append(f[1][-1])
    plt.plot(f[0], f[1], color="k", alpha=0.5)
tiempos = np.arange(0, end + 1)
plt.plot(tiempos, init*np.power(2, rate * tiempos), color="r")
plt.xlabel("tiempo")
plt.ylabel("tamaño de la población")
print "tamaño final promedio: ", np.mean(todoslostamaños)
plt.show()

respuesta a la excelente respuesta de whuber: no entiendo por qué tengo que especificar el tamaño de la población de antemano. mi simulación pretende preguntar: ¿en una cantidad determinada de tiempo, cuál será la variabilidad en el tamaño de la población asumiendo un crecimiento exponencial? (y no, ¿cuánto tiempo se necesita en promedio para hacer una población de tamaño $N$, que es lo que parece estar haciendo la simulación de whuber).

tampoco creo que sea correcto que "todavía no he hecho suficientes simulaciones para apreciar lo que te están diciendo". actualicé mi simulación para graficar 1000 corridas. como puedes ver, con mi condición de finalización basada en el tiempo, los resultados consistentemente subestiman el tamaño de la población de crecimiento exponencial determinista.

mi razonamiento era que si simulo $N$ corridas por una duración $t$, entonces a medida que $N \rightarrow \infty$, el tamaño promedio de la población que obtengo de mi simulación debería ser una estimación imparcial del tamaño de la población basada en el modelo de crecimiento exponencial determinista para $t$, es decir, $2^{rt}$. por ejemplo, si simulo 3 pasos de tiempo (comenzando con un solo individuo), esperaría que el tamaño de la población a veces sea mayor que 8 en la simulación y a veces menor que 8, y el promedio converja a 8. parece que esto debería ser cierto incluso si comienzo con una población muy pequeña, siempre y cuando simule suficientes corridas. ¿esto es incorrecto? ¿cuál es el error en este razonamiento? las simulaciones no respaldan esto a pesar de que esperaba que fuera cierto. parece que debe haber algún error en mi razonamiento y/o simulación.

actualización 2: corrección en la simulación donde el parámetro de escala de la distribución exponencial disminuye con el tamaño de la población (inversamente proporcional a él) y el tamaño de población inicial es 10. aún subestima enormemente el crecimiento exponencial.

Answer 1

El objetivo principal de la simulación es mostrarte que esa variación es realista. De hecho, no parece haber nada malo con tus resultados, excepto que aún no has hecho suficientes simulaciones para apreciar lo que te están diciendo. En este caso, los resultados son especialmente erráticos porque la población inicial es muy pequeña.

Corramos tu escenario 500 veces hasta una población de 300 (en lugar de media docena de veces hasta una población de 3):

Se ve más estable cuando comienzas con una población más grande:

Solo por diversión, aquí tienes una simulación similar para una población que crece desde un individuo hasta su capacidad de carga:

Usé R para estas simulaciones y gráficos, porque hace algo muy interesante. Dado que sabes de antemano que la población progresará en pasos enteros desde la población inicial hasta la final, puedes generar fácilmente la secuencia de valores de población antes de la simulación. Así, todo lo que queda por hacer es generar un conjunto de variaciones distribuidas exponencialmente con tasas determinadas por esa secuencia y acumularlos para simular los tiempos de nacimiento. R realiza esto con un solo comando (la línea que crea simulation a continuación). Toma alrededor de un segundo. Todo lo demás es solo especificación de parámetros y gráficos.

(Puedo usar un algoritmo tan simple porque estoy ejecutando estas simulaciones hasta un objetivo de población dado en lugar de hasta un punto final de tiempo dado. Obviamente, el modelo es el mismo; todo lo que difiere es cómo controlo la longitud de la simulación.)

rate <- 1
pop.0 <- 1
time.0 <- 0
k <-   0        # Capacidad de carga (usar 0 o negativo cuando no sea aplicable)
n.final <- 300  # ¡No debe exceder la capacidad!
#
# Precálculo: poblaciones y las tasas asociadas.
#
n <- pop.0:(n.final-1)
if (k <= 0) r <- rate * n else r <- rate * n * (1 - n/k)
#
# La simulación.
# Cada iteración se almacena como una columna del resultado.
#
simulation <- replicate(500, cumsum(c(time.0, rexp(length(n), r)))
#
# Graficar los resultados:
# Configurar, mostrar las curvas de crecimiento superpuestas, luego graficar una curva de referencia.
#
plot(range(simulation), c(pop.0, n.final), type="n", ylab="Población", xlab="Tiempo")
apply(simulation, 2, function(x) lines(x, c(n, n.final), col="#00000020"))
if (k <= 0) {curve(pop.0 * exp((x - time.0)*rate), add=TRUE, col="Red", lwd=2)} else
    curve(k*(1 - 1/(1+(pop.0/(k-pop.0))*exp(rate*(x-time.0)))), add=TRUE, col="Red", lwd=2)