3.1: Un modelo estocástico de crecimiento poblacional

El tamaño de la población\(N\) se considera ahora como una variable aleatoria discreta. Definimos la función de masa de probabilidad dependiente\(p_{N}(t)\) del tiempo de\(N\) ser la probabilidad de que la población sea de tamaño\(N\) en el tiempo\(t\). Ya que\(N\) debemos asumir uno de los valores de cero a infinito, tenemos

\[\sum_{N=0}^{\infty} p_{N}(t)=1 \nonumber \]

para todos\(t \geq 0 .\) De nuevo, deje\(b\) ser la tasa media de natalidad per cápita. Hacemos las aproximaciones simplificadoras de que todos los nacimientos son singletes, y que la probabilidad de que un individuo dé a luz es independiente de la historia pasada del parto. Entonces podemos interpretar\(b\) probabilísticamente suponiendo que como\(\Delta t \rightarrow 0\), la probabilidad de que un individuo dé a luz durante el tiempo\(\Delta t\) viene dada por\(b \Delta t\). Por ejemplo, si la tasa de natalidad per cápita promedio es de una descendencia cada 365 días del año, entonces la probabilidad de que un individuo dado dé a luz en un día determinado es\(1 / 365\). Como vamos a estar considerando el límite como\(\Delta t \rightarrow 0\), descuidamos las probabilidades de más de un nacimiento en la población en el intervalo de tiempo\(\Delta t\) ya que son de orden\((\Delta t)^{2}\) o superiores. Además, supondremos que at\(t=0\), se sabe que el tamaño de la población es\(N_{0}\), de modo que\(p_{N_{0}}(0)=1\), con todos los demás\(p_{N}^{\prime}\) s\(t=0\) iguales a cero.

Podemos determinar un sistema de ecuaciones diferenciales para la función de masa de probabilidad de la\(p_{N}(t)\) siguiente manera. Para que una población sea de tamaño\(N>0\) a la vez\(t+\Delta t\), o era de tamaño\(N-1\) a la vez\(t\) y se produjo un nacimiento, o bien era de tamaño\(N\) a la vez\(t\) y no hubo nacimientos; es decir

\[p_{N}(t+\Delta t)=p_{N-1}(t) b(N-1) \Delta t+p_{N}(t)(1-b N \Delta t) \nonumber \]

Restar\(p_{N}(t)\) de ambos lados, dividir por\(\Delta t\), y tomar los\(\Delta t \rightarrow 0\) resultados límite en las ecuaciones diferenciales de Kolmogorov hacia adelante,

\[\frac{d p_{N}}{d t}=b\left[(N-1) p_{N-1}-N p_{N}\right], \quad N=1,2, \ldots \nonumber \]

donde\(p_{0}(t)=p_{0}(0)\) ya que una población de tamaño cero permanece cero. Este sistema de ecuaciones diferenciales lineales acopladas de primer orden se puede resolver iterativamente. Primero revisamos cómo resolver una ecuación diferencial lineal de primer orden de la forma

\[\frac{d y}{d t}+a y=g(t), \quad y(0)=y_{0} \nonumber \]

donde\(y=y(t)\) y\(a\) es constante. Primero, buscamos un factor integrador\(\mu\) tal que

\[\frac{d}{d t}(\mu y)=\mu\left(\frac{d y}{d t}+a y\right) \nonumber \]

Diferenciar el lado izquierdo y multiplicar los resultados del lado derecho en

\[\frac{d \mu}{d t} y+\mu \frac{d y}{d t}=\mu \frac{d y}{d t}+a \mu y \nonumber \]

y términos de cancelación de rendimientos

\[\frac{d \mu}{d t}=a \mu . \nonumber \]

Podemos integrar esta ecuación con una condición inicial arbitraria, así que por simplicidad tomamos\(\mu(0)=1\). Por lo tanto,\(\mu(t)=e^{a t}\). Por lo tanto,

\[\frac{d}{d t}\left(e^{a t} y\right)=e^{a t} g(t) \nonumber \]

Integrando esta ecuación de 0 a t rendimientos

\[e^{a t} y(t)-y(0)=\int_{0}^{t} e^{a s} g(s) d s \nonumber \]

Por lo tanto, la solución es

\[y(t)=e^{-a t}\left(y(0)+\int_{0}^{t} e^{a s} g(s) d s\right) \nonumber \]

La ecuación diferencial de Kolmogorov hacia adelante (3.1.3) es de la forma (3.1.4) con\(a=b N\) y\(g(t)=b(N-1) p_{N-1}\). Con el tamaño de la población conocido por estar\(N_{0}\) en\(t=0\), las condiciones iniciales se pueden escribir sucintamente como\(p_{N}(0)=\delta_{N, N_{0}}\), donde\(\delta_{i j}\) está el delta de Kronecker, definido como

\[\delta_{i j}= \begin{cases}0, & \text { if } i \neq j \\[4pt] 1, & \text { if } i=j\end{cases} \nonumber \]

Por lo tanto, la integración formal de (3.1.3) usando (3.1.10) da como resultado

\[p_{N}(t)=e^{-b N t}\left[\delta_{N, N_{0}}+b(N-1) \int_{0}^{t} e^{b N s} p_{N-1}(s) d s\right] \nonumber \]

Las primeras soluciones de ahora se\((3.1.12)\) pueden obtener mediante integraciones sucesivas:

\[p_{N}(t)= \begin{cases}0, & \text { if } N<N_{0} \\[4pt] e^{-b N_{0} t}, & \text { if } N=N_{0} ; \\[4pt] N_{0} e^{-b N_{0} t}\left[1-e^{-b t}\right], & \text { if } N=N_{0}+1 ; \\[4pt] \frac{1}{2} N_{0}\left(N_{0}+1\right) e^{-b N_{0} t}\left[1-e^{-b t}\right]^{2}, & \text { if } N=N_{0}+2 ; \\[4pt] \cdots, & \text { if } \ldots\end{cases} \nonumber \]

A pesar de que no vamos a necesitar esto, para completitud doy la solución completa. Al definir el coeficiente binomial como el número de formas en las que se pueden seleccionar\(\mathrm{k}\) objetos de un conjunto de n objetos idénticos, donde el orden de selección es inmaterial, tenemos

\[\left(\begin{array}{l} n \\[4pt] k \end{array}\right)=\frac{n !}{k !(n-k) !} \nonumber \]

(leer como "\(n\)elegir\(k^{\prime \prime}\)). La solución general para\(p_{N}(t), N \geq N_{0}\), se sabe que es

\[p_{N}(t)=\left(\begin{array}{c} N-1 \\[4pt] N_{0}-1 \end{array}\right) e^{-b N_{0} t}\left[1-e^{-b t}\right]^{N-N_{0}} \nonumber \]

que los estadísticos llaman una distribución binomial negativa desplazada. La determinación de la evolución temporal de la función de masa probabilística de resuelve\(N\) completamente este problema estocástico.

De interés principal habitual es la media y varianza del tamaño de la población, y aunque ambos podrían calcularse en principio a partir de la función de masa de probabilidad, los calcularemos directamente a partir de la ecuación diferencial para\(p_{N}\). Las definiciones del tamaño medio de la población\(\langle N\rangle\) y su varianza\(\sigma^{2}\) son

\[\langle N\rangle=\sum_{N=0}^{\infty} N p_{N}, \quad \sigma^{2}=\sum_{N=0}^{\infty}(N-\langle N\rangle)^{2} p_{N} \nonumber \]

y haremos uso de la igualdad

\[\sigma^{2}=\left\langle N^{2}\right\rangle-\langle N\rangle^{2} \nonumber \]

Multiplicando la ecuación diferencial (3.1.3) por la constante\(N\), sumando\(N\) y usando\(p_{N}=0\) for\(N<N_{0}\), obtenemos

\[\frac{d\langle N\rangle}{d t}=b\left[\sum_{N=N_{0}+1}^{\infty} N(N-1) p_{N-1}-\sum_{N=N_{0}}^{\infty} N^{2} p_{N}\right] \nonumber \]

Ahora, escribe\(N(N-1)=(N-1)(N-1+1)=(N-1)^{2}+(N-1)\), para que el primer término en el lado derecho sea

\[\begin{aligned} \sum_{N=N_{0}+1}^{\infty} N(N-1) p_{N-1} &=\sum_{N=N_{0}+1}^{\infty}(N-1)^{2} p_{N-1}+\sum_{N=N_{0}+1}^{\infty}(N-1) p_{N-1} \\[4pt] &=\sum_{N=N_{0}}^{\infty} N^{2} p_{N}+\sum_{N=N_{0}}^{\infty} N p_{N} \end{aligned} \nonumber \]

donde la segunda igualdad se obtuvo desplazando el índice de suma a la baja en uno. Por lo tanto, encontramos la ecuación familiar del crecimiento maltusiano

\[\frac{d\langle N\rangle}{d t}=b\langle N\rangle . \nonumber \]

Junto con la condición inicial\(\langle N\rangle(0)=N_{0}\), podemos encontrar la solución

\[\langle N\rangle(t)=N_{0} e^{b t} \nonumber \]

Se procede de manera similar a encontrar\(\sigma^{2}\) determinando primero la ecuación diferencial para\(\left\langle N^{2}\right\rangle\). Multiplicar la ecuación diferencial para\(p_{N},(3.1.3)\), por\(N^{2}\) y sumando sobre\(N\) los resultados en

\[\frac{d\left\langle N^{2}\right\rangle}{d t}=b\left[\sum_{N=N_{0}+1}^{\infty} N^{2}(N-1) p_{N-1}-\sum_{N=N_{0}}^{\infty} N^{3} p_{N}\right] \nonumber \]

Aquí, escribimos\(N^{2}(N-1)=(N-1)(N-1+1)^{2}=(N-1)^{3}+2(N-1)^{2}+(N-1)\). Procediendo de la misma manera que anteriormente desplazando el índice a la baja, obtenemos

\[\frac{d\left\langle N^{2}\right\rangle}{d t}-2 b\left\langle N^{2}\right\rangle=b\langle N\rangle \nonumber \]

Dado que\(\langle N\rangle\) se conoce, (3.1.22) es una ecuación de primer orden, lineal, no homogénea para\(\left\langle N^{2}\right\rangle\), la cual puede resolverse usando un factor integrador. La solución obtenida usando (3.1.10) es

\[\left\langle N^{2}\right\rangle=e^{2 b t}\left(N_{0}^{2}+b \int_{0}^{t} e^{-2 b s}\langle N\rangle(s) d s\right) \nonumber \]

con\(\langle N\rangle\) dado por (3.1.20). Al realizar la integración, obtenemos

\[\left\langle N^{2}\right\rangle=e^{2 b t}\left[N_{0}^{2}+N_{0}\left(1-e^{-b t}\right)\right] \nonumber \]

Finalmente, usando\(\sigma^{2}=\left\langle N^{2}\right\rangle-\langle N\rangle^{2}\), obtenemos la varianza. Así llegamos a nuestro resultado final para la media poblacional y varianza:

\[\langle N\rangle=N_{0} e^{b t}, \quad \sigma^{2}=N_{0} e^{2 b t}\left(1-e^{-b t}\right) \nonumber \]

El coeficiente de variación\(c_{v}\) mide la desviación estándar relativa a la media, y aquí viene dada por

\[\begin{aligned} c_{v} &=\sigma /\langle N\rangle \\[4pt] &=\sqrt{\frac{1-e^{-b t}}{N_{0}}} . \end{aligned} \nonumber \]

Para grandes\(t\), el coeficiente de variación por lo tanto va como\(1 / \sqrt{N_{0}}\), y es pequeño cuando\(N_{0}\) es grande. En la siguiente sección, determinaremos la forma limitante de la distribución de probabilidad para grandes\(N_{0}\), recuperando tanto el modelo determinista como una aproximación del modelo gaussiano.