Modelación estocástica
4 octubre 2021
Abstract
La incertidumbre juega un papel cada vez mayor en la forma en que se entiende al mundo, ya sea el clima del planeta, la propagación de una enfermedad infecciosa o los resultados de la última encuesta de noticias. Un proceso estocástico, también llamado proceso aleatorio, es simplemente uno en el que los resultados son inciertos. Por el contrario, en un sistema determinista no hay aleatoriedad. En un sistema determinista, la misma salida siempre se produce a partir de un entrada dada. Las funciones y las ecuaciones diferenciales se utilizan normalmente para describir procesos deterministas. Las variables aleatorias y las distribuciones de probabilidad son las componentes básicas de sistemas estocásticos.
1 Modelos deterministas y estocásticos
Considérese que la cantidad de bacterias que crecen en un cultivo hasta que se agota su fuente de alimento exhibe un crecimiento exponencial. Un modelo de crecimiento determinista común es afirmar que la población de bacterias crece a tasa fija, digamos \(20%\) por minuto. El modelo determinista no aborda la incertidumbre presente en la tasa de reproducción de organismos individuales. Tal incertidumbre se puede capturar mediante el uso de un marco estocástico en el que los tiempos hasta que las bacterias se reproducen se modelan mediante variables aleatorias.
Un modelo de crecimiento estocástico simple es asumir que los tiempos hasta que las bacterias individuales se reproducen son variables aleatorias exponenciales independientes, en este caso con un parámetro de tasa de \(0.20\).
En muchos procesos biológicos, la distribución exponencial es una opción común para modelar los tiempos de nacimientos y muertes. En el modelo determinista, cuando el tamaño de la población es \(n\), el número de bacterias aumenta en (\(0.20\)) \(n\) en \(14\) minutos.
De manera similar, para el modelo estocástico, después de que surgen \(n\) bacterias, el tiempo hasta que la siguiente bacteria se reproduce tiene una distribución de probabilidad exponencial con una tasa (\(0.20\)) \(n\) por minuto. Si bien el resultado de un sistema determinista es fijo, el resultado de un proceso estocástico es incierto.
La dinámica de un proceso estocástico se describe mediante variables aleatorias y distribuciones de probabilidad. En el modelo de crecimiento determinista, se puede decir con certeza cuántas bacterias están presentes después de \(t\) minutos. Para el modelo estocástico, las preguntas de interés pueden incluir:
- ¿Cuál es el número promedio de bacterias presentes en el momento \(t\)?
- ¿Cuál es la probabilidad de que el número de bacterias supere algún umbral después de \(t\) minutos?
- ¿Cuál es la distribución del tiempo que tarda la cantidad de bacterias en duplicar su tamaño?
El crecimiento exponencial se puede modelar mediante la ecuación \(y (t) = ae^{bt}, a, b> 0, t \in [0,\infty)\) y se comparan con el modelo de crecimiento exponencial, \(y (t) = e^t\). A continuación se emplearán rutinas en Octave basadas en Allen (2010) y para R en Dobrow (2016).
Una parte importante del modelado es la simulación numérica. Se pueden utilizar muchos lenguajes de programación diferentes para simular la dinámica de un modelo estocástico. La salida de un programa en Octave (Eaton et al. 2020) para tres realizaciones estocásticas del proceso de nacimiento simple cuando \(b = 1\), \(a = 1\) y \(p_1 (0) = 1\) se grafican en la Figura 1. También se representa gráficamente el modelo de crecimiento exponencial determinista correspondiente, \(n (t) = e^t\).
Figure 1: Realizaciones estocásticas del proceso de nacimiento simple cuando b = 1, a = 1
El modelo de crecimiento exponencial determinista, n(t) = e^t, es la curva discontinua.
Las nociones computacionales de una función de distribución y las densidades se dan en muchas fuentes, por ejemplo en R (R Core Team 2021) y RStudio (RStudio Team 2021), se cuenta con Dutang and Kiener (2018), Robert A. Rigby et al. (2020), R. A. Rigby and Stasinopoulos (2005), Dutang, Goulet, and Pigeon (2008) y Rassel and Crépeault-Cauchon (2020).
Para la mayoría de las distribuciones clásicas, R base proporciona funciones de distribución de probabilidad (p), funciones de densidad (d), funciones de cuantiles (q) y generación de números aleatorios (r).
Como ejemplo, en Octave se pueden generar gráficas para la poisson, como se muestra en en la Figura 2.
Figure 2: Función de masa de probabilidad Poisson con parámetro λ = 3
Los modelos estocásticos para la propagación de enfermedades infecciosas y el desarrollo de epidemias son relevantes debido a la aleatoriedad inherente a los contactos de persona a persona y las fluctuaciones de la población. El modelo SIR (Susceptible-Infected-Removed) es un marco básico que se ha aplicado a la propagación del sarampión y otras enfermedades infantiles. En el tiempo \(t\), sea \(S_t\) el número de personas susceptibles a una enfermedad, \(I_t\) el número de infectados y \(R_t\) el número recuperado y, en adelante, inmunes a la infección. Los individuos de la población pasan de ser susceptibles a posiblemente infectados a recuperarse (S → I → R).
El modelo SIR determinista se deriva de un sistema de tres ecuaciones diferenciales no lineales, que modelan las interacciones y la tasa de cambio en cada subgrupo.
Un modelo SIR estocástico en tiempo discreto fue introducido en la década de \(1920\) por Lowell Reed y Wade Frost de la Universidad Johns Hopkins. En el modelo de Reed-Frost, cuando un individuo susceptible entra en contacto con alguien que está infectado, existe una probabilidad fija \(z\) de que se infecte.
Suponiendo que cada persona susceptible está en contacto con todas las que están infectadas. Sea \(p\) la probabilidad de que un individuo susceptible esté infectado en el tiempo \(t\). Esto es igual a \(1\) menos la probabilidad de que la persona no esté infectada en el momento \(t\), lo que ocurre si no está infectada por ninguna de las personas ya infectadas, de las cuales hay \(I_t\).
Esto da
\[\begin{equation} p = 1 − (1 − z)^{I_t}. \tag{1} \end{equation}\]
La evolución de la enfermedad se modela en tiempo discreto, donde una unidad de tiempo es el período de incubación, también el tiempo de recuperación, de la enfermedad.
El modelo se puede describir con una analogía del lanzamiento de una moneda. Para encontrar \(I_{t + 1}\), el número de individuos infectados en el tiempo \(t + 1\), se lanzan \(S_t\) monedas (una por cada susceptible), donde la probabilidad de que salga cara por cada moneda es la probabilidad de infección \(p\). Entonces, el número de individuos recién infectados es el número de monedas que caen cara. El número de caras en \(n\) lanzamientos de monedas independientes con probabilidad de caras \(p\) tiene una distribución binomial con parámetros \(n\) y \(p\). En otras palabras, \(I_{t + 1}\) tiene una distribución binomial con \(n = S_t\) y \(p = 1 - (1 - z)^{I_t}\).
Habiendo encontrado el número de individuos infectados en el tiempo \(t + 1\), el número de personas susceptibles disminuye por el número de infectados. Eso es,
\[\begin{equation} S _{t+1} = S _t − I _t+1. \tag{2} \end{equation}\]
Aunque el modelo Reed-Frost no es fácil de analizar con exactitud, es sencillo de simular. Las gráficas de la Figura 3 se obtuvieron simulando el proceso asumiendo una población inicial de \(3\) infectados y \(400\) individuos susceptibles, con una probabilidad de infección individual \(z = 0.004\). El número de infectados se representa en \(20\) unidades de tiempo.
# Simulación ReedFrost
par(mfrow=c(2,2))
set.seed(12)
for (k in 1:4) {
z = 0.004
steps <- 21
sus <- numeric(steps)
inf <- numeric(steps)
sus[1] <- 400
inf[1] <- 3
for (t in 1:(steps-1)) {
inf[t+1] <- rbinom(1,sus[t],1-(1-z)^inf[t])
sus[t+1] <- sus[t] - inf[t+1]
}
plot((0:(steps-1)),inf,type="o",yaxt="n",ylim=c(0,36),xlab="Time",ylab="Infected")
axis(2,at=c(0,10,20,30))
}
Figure 3: Cuatro realizaciones del modelo Reed-Frost
#seed: 122 Modelación con bases de datos
Además de las funcionalidades básicas de R, se cuenta con la manipulación del tidyverse (Wickham et al. 2019).
En el proceso de análisis de epidemias desde el AED hasta la difusión mediante un estudio de caso de sarampión, la librería EpiCompare (Gallagher and LeRoy 2021) proporciona herramientas fáciles de usar para fomentar la comparación y evaluación de epidemias y modelación epidemiológica. Todas las herramientas intentan adherirse al estilo de tidyverse / ggplot2 (Wickham 2016) para mejorar la facilidad de uso.
2.1 Pipeline del análisis de epidemias
Figure 4: Flujo de trabajo en EpiCompare (Gallagher and LeRoy, 2021)
2.2 Pre-procesamiento y AED
El conjunto de datos Hagelloch sigue el curso de una epidemia de sarampión en Hagelloch, Alemania, desde el \(30\) de octubre de \(1861\) (día \(0\)) hasta el \(24\) de enero de \(1862\), cubriendo un período de \(87\) días.
Junto con las paperas, la rubéola y la varicela, el sarampión es una enfermedad infantil altamente infecciosa. El sarampión se transmite de persona a persona a través del aire contaminado o una superficie infectada, es quizás la enfermedad más contagiosa del planeta, con un número de reproducción estimado de alrededor de \(R0 = 19\), lo que significa que cuando una persona infecciosa se introduce en una población totalmente susceptible, infectará a otras \(19\) en promedio.
Se pueden ver los datos sin procesar utilizando data(hagelloch_raw) y ver los detalles con ?hagelloch_raw. Una parte de la base de datos brutos se muestra en la Tabla 1.
library(tidyverse)
library(EpiCompare)
library(knitr)
library(kableExtra)
library(RColorBrewer)
hagelloch_raw %>% select(PN, NAME, AGE, SEX, HN, PRO, ERU, IFTO) %>%
head(10) %>% kable( caption = "Hagelloch data sin procesar.") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| PN | NAME | AGE | SEX | HN | PRO | ERU | IFTO |
|---|---|---|---|---|---|---|---|
| 1 | Mueller | 7 | female | 61 | 1861-11-21 | 1861-11-25 | 45 |
| 2 | Mueller | 6 | female | 61 | 1861-11-23 | 1861-11-27 | 45 |
| 3 | Mueller | 4 | female | 61 | 1861-11-28 | 1861-12-02 | 172 |
| 4 | Seibold | 13 | male | 62 | 1861-11-27 | 1861-11-28 | 180 |
| 5 | Motzer | 8 | female | 63 | 1861-11-22 | 1861-11-27 | 45 |
| 6 | Motzer | 12 | male | 63 | 1861-11-26 | 1861-11-29 | 180 |
| 7 | Schneck | 6 | male | 23 | 1861-11-24 | 1861-11-28 | 42 |
| 8 | Schneck | 10 | male | 69 | 1861-11-21 | 1861-11-26 | 45 |
| 9 | Schneck | 13 | male | 69 | 1861-11-26 | 1861-11-30 | 182 |
| 10 | Schneck | 7 | female | 31 | 1861-11-21 | 1861-11-25 | 45 |
donde PRO es la fecha en la que comenzaron los pródromos (es decir, los síntomas iniciales) y ERU es la fecha de la erupción del sarampión. Estos datos también incluyen la variable IFTO, que es la identificación del infectante (sospechosa), lo que hace que estos datos sean muy informativos e inusuales. Solo el \(2%\) de la población infantil infectada de un total de \(188\) no tienen una identificación de infectante reportada.
La base también incluyen tI y tR que son tiempos estimados de infección y recuperación. La duración de la infección infantil de la Figura 5, colorea líneas por hogar (con colores reciclados) y organizando a infantes por hogar. A partir de esto, se ve que la infección parece transmitirse a través del hogar en un período corto de tiempo donde casi siempre hay superposición dentro del hogar. También, al colorear las líneas por grado escolar se muestra que la primera clase parece haberse infectado aproximadamente al mismo tiempo.
hagelloch_raw %>% arrange(HN, NAME) %>%
mutate(new_id = paste(HN, NAME, PN, sep = "-")) %>%
group_by(HN) %>%
mutate(new_id = forcats::fct_reorder(new_id, -tI)) %>%
ggplot2::ggplot(aes(y = new_id, col = factor(HN))) +
scale_colour_manual(values=rep(brewer.pal(5,"Set1"),
times= ceiling(length(unique(hagelloch_raw$HN)) / 5)),
guide = FALSE) +
geom_errorbarh(aes(xmax = tR, xmin = tI), size = 1) + theme_bw() +
labs(x = "Date",
y = "Child ID",
title = "Infection time of children",
subtitle = "Colored by household ID (with recycled colors)") +
theme(axis.text.y = element_text(size = 4))
Figure 5: Duración de infección infantil por sarampión
Y por grado escolar en la Figura 6.
hagelloch_raw %>% arrange(CL) %>%
mutate(new_id = paste(CL, PN, sep = "-")) %>%
group_by(CL) %>%
mutate(new_id = forcats::fct_reorder(new_id, -tI)) %>%
ggplot2::ggplot(aes(y = new_id, col = CL))+
geom_errorbarh(aes(xmax = tR, xmin = tI), size = 1) + theme_bw() +
labs(x = "Date",
y = "Child ID",
title = "Infection time of children",
color = "School Grade") +
viridis::scale_color_viridis(discrete = TRUE, end = .8) +
theme(axis.text.y = element_text(size = 4))
Figure 6: Duración de infección infantil por sarampión
Información por grado escolar.
Al contar las infecciones producidas, mostradas en la Figura 7, por lo que se sospecha que una infeccción se esparcec a \(30\) demás susceptibles, la mayoría de las cuales son otras personas en su misma clase.
hagelloch_raw$inf_trans <- sapply(hagelloch_raw$PN, function(x){
sum(hagelloch_raw$IFTO == x, na.rm = TRUE)}
)
hagelloch_raw %>% filter(inf_trans > 0) %>%
ggplot() + geom_col(aes(x = reorder(PN, inf_trans), y = inf_trans,
fill = CL)) +
coord_flip() +
viridis::scale_fill_viridis(discrete = TRUE, end = .8) +
labs(x = "ID", y = "Number of (suspected) transmissions",
title = "(Suspected) Number of infections caused by child") +
theme_bw()
Figure 7: Duración de infección infantil por sarampión
Información por grado escolar.
Ver información a nivel agente es interesante, pero se puede visualizar información a nivel general, es decir, el número total de infantes en un estado de infección determinado a lo largo del tiempo. Con la función EpiCompare::agent_to_aggregate() se convierten los datos del nivel del agente en una vista agregada en la Tabla 2. La Figura 8 presenta el número de infantes infecciosos a lo largo del tiempo.
aggregate_hag <- hagelloch_raw %>%
agents_to_aggregate(states = c(tI, tR)) %>%
rename(time = t, S = X0, I = X1, R = X2)
aggregate_hag %>% head %>%
kable( caption = "Hagelloch data vista agregada.") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| time | S | I | R |
|---|---|---|---|
| 0 | 187 | 1 | 0 |
| 1 | 187 | 1 | 0 |
| 2 | 186 | 2 | 0 |
| 3 | 186 | 2 | 0 |
| 4 | 186 | 2 | 0 |
| 5 | 186 | 2 | 0 |
aggregate_hag %>%
ggplot(aes(x = time, y = I)) +
geom_line(col = "red") +
geom_point(col = "red") +
theme_bw() +
labs( x= "Time after initial infection (days)",
y = "# of infected children",
title = "Number of infectious children vs. time",
subtitle = "Hagelloch, Germany 1861-1862")
Figure 8: Infantes infecciosos a lo largo del tiempo
EpiCompare::agent_to_aggregate() también funciona en data frames agrupados, lo que significa que se pueden agregar infantes por grupos como su clase, en la Tabla 3. La Figura 9 presenta el número de infantes infecciosos a lo largo del tiempo por clase.
aggregate_hag_class <- hagelloch_raw %>%
group_by(CL) %>%
agents_to_aggregate(states = c(tI, tR)) %>%
rename(time = t, S = X0, I = X1, R = X2)
aggregate_hag_class %>% head %>%
kable( caption = "Hagelloch data vista agregada. Agrupación por clase.") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| CL | time | S | I | R |
|---|---|---|---|---|
| 1st class | 0 | 30 | 0 | 0 |
| 1st class | 1 | 30 | 0 | 0 |
| 1st class | 2 | 30 | 0 | 0 |
| 1st class | 3 | 30 | 0 | 0 |
| 1st class | 4 | 30 | 0 | 0 |
| 1st class | 5 | 30 | 0 | 0 |
aggregate_hag_class %>%
ggplot(aes(x = time, y = I, group = CL,
col = CL)) +
geom_line() +
geom_point() +
theme_bw() +
viridis::scale_color_viridis(discrete = TRUE, end = .8) +
labs( x= "Time after initial infection (days)",
y = "# of infected children",
title = "Number of infectious children vs. time",
subtitle = "Hagelloch, Germany 1861-1862")
Figure 9: Infantes infecciosos a lo largo del tiempo
Por clase.
Lo que sugiere que la primera clase pudo haberse infectado antes que las otras dos.
En resumen, se presentó la base de datos de sarampión de Hagelloch a nivel de agente y a nivel agregado. La función EpiCompare::agents_to_aggregate() preprocesa datos y a un formato para ver fácilmente información de nivel agregado que sea compatible con ggplot2 y otras funciones de estilo tidyverse.
2.3 Modelos
El modelo SIR se basa en una serie de ecuaciones diferenciales ordinarias (EDO). La idea es que los individuos fluyan de un estado (o compartimento) discreto al siguiente. Estas transiciones de compartimento a compartimento están definidas por EDO:
\[\begin{equation} \begin{gathered} S'(t)=−βS(t)I(t) \\ I'(t)=βS(t)I(t)−γI(t) \\ R'(t)=γI(t), \end{gathered} \tag{3} \end{equation}\]
con restricciones \(S'(t) + I'(t) + R'(t) ≡0\) y condiciones inicicales \((S(0),I(0),R(0))=(s(0),i(0),r(0))\).
En este conjunto de EDO, la tasa instantánea de cambio en el número de susceptibles es igual al parámetro \(β\) veces el número de susceptibles y el número de infecciosos. Otra forma de pensar en esto es que los susceptibles e infecciosos interactúan de manera homogénea entre sí y se infectan a un ritmo \(β\). La tasa instantánea de cambio del número de recuperados es igual al parámetro \(γ\) veces el número de infecciosos, o que los infecciosos se están recuperando a un ritmo \(γ\). Entonces finalmente, \(I′( t ) = -S′( t ) -R′( t ),\) ya que la suma de las derivadas es cero.
Un modelo SIR de tiempo discreto y determinista discretiza los pasos de tiempo. Asumiendo el paso de tiempo \(Δt = 1\), y las ecuaciones en diferencias para \(t=1,...,T\) son:
\[\begin{equation} \begin{gathered} -βS(t-1)I(t-1)=[S(t)-S(t-1)] \\ [I(t)-I(t-1)]=βS(t-1)I(t-1)-γI(t-1)\\ [S(t)-S(t-1)]=γI(t-1), \end{gathered} \tag{4} \end{equation}\]
que son las aproximaciones de primer orden de las EDO.
Hacer una función objetiva para minimizar el modelo SIR determinista y de tiempo continuo da mejor ajuste al considerar la suma del error al cuadrado sobre los datos observados y las estimaciones:
\[\begin{equation} \begin{gathered} (\hat{β},\hat{y}) & = arg \min_{β, γ}∑_{t = 0}^T ( S( t ) - s ( t ))^2 + ( I( t ) - i ( t ))^2+ \\ & ( R ( t ) - r ( t ))^2, \end{gathered} \tag{5} \end{equation}\]
donde \(S(t),I(t),\) y \(R(t)\) son estimaciones producidas por el modelo y \(s(t),i(t),r(t)\) es información observable.
Transformamos este modelo SIR determinista de tiempo continuo en uno estocástico de tiempo discreto. El primer paso es cambiar el \(Ds\) para \(Δs\), o más formalmente usando una aproximación de primer orden de las EDO.
El primer paso para transformar eL modelo SIR determinista de tiempo continuo en uno estocástico de tiempo discreto es cambiar el \(ds\) por un \(Δs\), o más formalmente hacer uso de una aproximación de primer orden de las EDO.
\[\begin{equation} \begin{gathered} ΔS(t)Δt=S(t)−S(t−1)Δt=−βS(t−1)I(t−1) \\ ΔI(t)Δt=I(t)−I(t−1)Δt=βS(t−1)I(t−1)−γI(t−1) \\ ΔR(t)Δt=R(t)−R(t−1)Δt=γI(t−1). \end{gathered} \tag{6} \end{equation}\]
Asumiendo \(Δt = 1\), el siguiente paso es introducir un error aleatorio en el modelo. Específicamente, el error de proceso, que es un error inherente al proceso de transmisión, en contraposición al error introducido por el error sistémico en la medición de los datos.
Hay varias formas de modelar el número de individuos que se mueven de un estado a otro de un paso de tiempo al siguiente, pero algunas de las más comunes son las variables aleatorias binomiales y poisson. Considerando el proceso con variables aleatorias binomiales, condicionadas a los valores del estado anterior.
Dados los valores iniciales \(X_S(0),X_I(0),X_R(0)\), el número de personas que se mueven de un estado a otro de un tiempo \(t - 1\) para \(t\) con \(t > 0\) vienen dadas por las siguientes ecuaciones:
\[\begin{equation} \begin{gathered} Z_{SI}(t)|X_S(t−1),X_I(t−1)∼Binomial(X_S(t−1),βX_I(t−1)) \\ Z_{IR}(t)|X_I(t−1)∼Binomial(X_I(t−1),γ). \end{gathered} \tag{7} \end{equation}\]
El número de personas que pasan del estado susceptible al estado infeccioso desde el momento \(t - 1\) para \(t\), \(Z_{SI}( t )\), está eligiendo entre \(X_S( t - 1 )\) susceptibles y con probabilidad de \(βI( t - 1 )\) toma a varios de ellos para que se vuelvan infecciosos. Similarmente ,\(Z_{IR}( t )\) elige de \(X_I( t - 1 )\) individuos infectados con probabilidad \(γ\) para recuperarse. Esto significa que \(βI( t - 1 )\) y \(γ\) tienen que estar restringidas entre cero y uno, para tener una distribución correctamente especificada.
El nuevo número en cada estado viene dado por lo siguiente:
\[\begin{equation} \begin{gathered} X_S(t)=X_S(t−1)−Z_{SI} \\ X_I(t)=X_I(t−1)+Z_{SI}−Z_{IR} \\ X_R(t)=X_R(t−1)+Z_{IR}. \end{gathered} \tag{8} \end{equation}\]
En este enfoque \(E[ (X_S( t ) ,X_I( t ) ,X_R( t ) ) ] = ( S( t ) , I( t ) , R ( t ) )\) para todo \(t\), y se puede calcular la varianza.
Denotando al conjunto de datos observados (el número de individuos en cada estado en el momento \(t\), \(t = 0 ,..., T\)) como \(\mathbf{x}\) y a la probabilidad de infectarse de vez \(t - 1\) para \(t\) como \(p_t= βX_I( t )\)., entonces la verosimilitud es entonces proporcional a:
\[\begin{equation} \begin{gathered} \mathcal{L} (β, γ; \mathbf{x} ) & = ∏_{t = 1}^T p^{S( t - 1 ) - S( t )}_{t - 1}[ 1 -p_{t - 1}]^{S( t )} \times \\ & γ^{R ( t ) - R ( t - 1 )}( 1 - γ)^{I( t - 1 ) - ( R ( t ) - R ( t - 1 ) )}, \end{gathered} \tag{9} \end{equation}\]
donde los exponentes son cantidades no aleatorias, sino que están completamente determinadas por los valores iniciales y los de \(β\) y \(γ\). La log-verosimilitud es entonces:
\[\begin{equation} \begin{gathered} \mathcal{l} ( β, γ, x ) & = ∑_{t = 1}^T( S( t - 1 ) - S( t ) ) \log (p_{t - 1}) + \\ & S( t ) \log ( 1 -p_{t - 1}) + ( R ( t ) - R ( t - 1 ) ) \log ( γ) + \\ & ( I( t - 1 ) - ( R ( t ) - R ( t - 1 ) ) ) \log ( 1 - γ). \end{gathered} \tag{10} \end{equation}\]
Una aproximación a la probabilidad logarítmica anterior se obtiene sustituyendo los datos \(x_S( t ), x_I( t ),\) y \(x_R( t )\) por \(S( t ) , I( t ), y R ( t )\), respectivamente.
Así, se toman los datos SIR agregados hagelloch_raw y se define un modelo SIR determinista y uno de tiempo continuo para encontrar las mejores estimaciones en la Tabla 4.
# datos agregados
aggregate_hag <- hagelloch_raw %>%
agents_to_aggregate(states = c(tI, tR),
min_max_time = c(0, 55)) %>%
rename(time = t, S = X0, I = X1, R = X2)
# calcular SIR log verosimilitud de eq. 10
# @param par: parametros a optimizar
# @param data: un data frame con columnas tiempo, S, I, y R
sir_loglike <- function(par, data){
data$beta <- par[1]
data$gamma <- par[2]
data$N <- data$S + data$I + data$R
data_new <- data %>%
dplyr::mutate(delta_S = c(NA, diff(S)),
delta_R = c(NA, diff(R)),
prev_S = dplyr::lag(S),
prev_I = dplyr::lag(I),
prev_R = dplyr::lag(R),
prob_inf = beta / N * prev_I,
loglike_SI_init = (-delta_S) * log(prob_inf) +
S * log(1 - prob_inf),
loglike_SI = ifelse(prev_I == 0, 0,
loglike_SI_init),
loglike_IR = (delta_R) * log(gamma) + (prev_I - delta_R) * log(1 - gamma),
loglike = loglike_SI + loglike_IR
) %>%
dplyr::filter(time != 0)
return(sum(data_new$loglike))
}
init_params <- c(.1, .05)
best_params <- optim(par = init_params, fn = sir_loglike,
control = list(fnscale = -1), # switches from minimizing function to maximizing it
data = aggregate_hag,
hessian = TRUE)
#print(best_params$par, digits = 2)
best_params$par %>% kable( caption = "Mejores estimaciones." , col.names = "best_params$par") %>%
kable_styling(bootstrap_options = c("striped", "hover"))| best_params$par |
|---|
| 0.362442 |
| 0.126167 |
Con la función simulate_SIR_agents() se simulan datos SIR a nivel individual (o agente), los agrega en formato SIR y los traza en la vista estándar de porcentaje en estado frente a tiempo (Figuras 10 y 11).
n_sims <- 1
n_time_steps <- 100
beta <- .1
gamma <- .03
init_SIR <- c(950, 50, 0)
out <- EpiCompare::simulate_SIR_agents(n_sims = n_sims,
n_time_steps = n_time_steps,
beta = beta, gamma = gamma,
init_SIR = init_SIR)
df_ave <- out %>% agents_to_aggregate(states = c(tI, tR)) %>%
rename(S = "X0", I = "X1", R = "X2") %>%
group_by(t) %>%
summarize(S = mean(S), I = mean(I), R = mean(R),
.groups = "drop")
df_ave_lines <- df_ave %>% pivot_longer(names_to = "State",
values_to = "Value",
cols = c(S, I, R))
ggplot(data = df_ave_lines,
aes(x = t, y = Value, group = State, col = State)) +
geom_line()
Figure 10: Simulación de valores I, R y S a lo largo del tiempo
ggplot ( data = df_ave_lines , aes ( x = t , y = Value )) +
facet_grid ( State ~ . ) +
geom_line ()
Figure 11: Simulación de valores I, R y S a lo largo del tiempo
Vista en facet.
Y se visualiza esta simulación con una sola pipeline en una gráfica ternaria (Figura 12).
ggplot ( data = df_ave , aes (x = S , y = I , z = R )) +
geom_path () + coord_tern ()
Figure 12: Simulación de valores I, R y S a lo largo del tiempo
Vista ternaria.
3 Simulación
La simulación es una herramienta eficaz para analizar modelos de probabilidad, así como para facilitar la comprensión de conceptos en probabilidad y estadística.
El paquete Symbulate (Ross and Sun 2019; Ross 2020, https://github.com/dlsun/symbulate) de Python (Van Rossum and Drake 2009) proporciona una herramienta fácil de usar para realizar simulaciones y es totalmente compatible con Jupyter (Kluyver et al. 2016), que también proporciona una interfaz fácil de usar que admite programación interactiva y documentación reproducible.
La sintáxis refleja el lenguaje de probabilidad en el sentido de que los objetos principales en la simulación sean los mismos que los componentes principales de un modelo de probabilidad:
Un espacio de probabilidad que consta del espacio muestral (Ω) de posibles resultados (ω∈Ω) y una medida de probabilidad (P),
Variables aleatorias (funciones X: Ω↦R) o procesos estocásticos (funciones X: T × Ω↦R)
Para un conjunto de índices (de tiempo) T definido en el espacio de probabilidad, posiblemente a través de transformaciones de variables aleatorias o procesos, y eventos relacionados (subconjuntos A⊂Ω para los cuales se define P (A)), y en particular, condicionamiento en eventos.
Con las componentes especificadas, se simula muchas veces a partir del modelo de probabilidad y se resumen esos resultados.
3.1 Espacios y variables aleatorias
El espacio de probabilidad en muchas situaciones elementales se puede definir mediante un modelo de caja. Como ejemplo al evaluar la intensidad de cierta infección al seleccionar a dos personas susceptibles, sea \(X\) la suma de estas dos variables aleatorias vistas como el lanzamiento de una moneda (\(1\)=Infección, \(0\)=No Infección), y sea \(Y\) la mayor de las variables (vista como el mayor de dos lanzamientos o el valor común en empate). Al definir un espacio de probabilidad \(P\) se pueden simular las \(4\) parejas ordenadas igualmente probables.
from symbulate import *
from matplotlib import *
P = BoxModel([0, 1], size = 2, replace = True)
P.sim(4)Con la especificación dada, una variable aleatoria es una función que asigna los resultados de un espacio de probabilidad a números reales: \(X(ω)\) para \(ω ∈ Ω\), el espacio muestral. Una VA en Symbulate se especifica mediante el espacio de probabilidad en el que se define y la función de mapeo.
X = RV(P, sum)
Y = RV(P, max) Al definir y manipular espacios de probabilidad y variables aleatorias como objetos abstractos, independientemente de las realizaciones simuladas, se facilita la comprensión de estos conceptos distintos.
Dado que una variable aleatoria \(X\) es una función, cualquier VA se puede llamar como una función para devolver su valor \(X(ω)\) para un resultado particular \(ω\) en el espacio de probabilidad.
outcome = (0, 1) # una pareja inicial
X(outcome), Y(outcome) Luego, se simulan \(10000\) valores de la variable aleatoria \(Y\), se almacenan los resultados en una variable \(y\) y se resumen los valores obtenidos y sus frecuencias relativas, reflejando la notación estándar y enfatizando la distinción explícita entre la variable aleatoria \(Y\) (letra mayúscula) y los valores realizados de ella \(y\) (letra minúscula).
y = Y.sim(10000)
y.tabulate(normalize = True)
y.plot()
y.count_leq(3)/10000, y.mean(), y.sd()Los valores de múltiples variables aleatorias se pueden simular simultáneamente, es decir, Para los mismos resultados simulados en el espacio de probabilidad. El siguiente código produce una gráfica que ilustra las distribuciones conjunta y marginal aproximadas de X e Y obtenidas por simulación.
Para los mismos resultados simulados en el espacio de probabilidad, los valores de múltiples variables aleatorias se pueden simular simultáneamente, como las distribuciones conjunta y marginal aproximadas de \(X\) y \(Y\) obtenidas por simulación.
3.2 Independencia
En Symbulate, la independencia está representada por el asterisco *, que refleja que bajo independencia, los objetos conjuntos (por ejemplo, probabilidad, cdf, pdf) son productos de los objetos marginales correspondientes.
El espacio de probabilidad \(P\) corresponde a un solo lanzamiento de una moneda de dos lados justa, y suponiendo que ahora hay dos cepa de la infección original, sea \(Q\) el espacio de un solo lanzamiento de un dado de cuatro lados justo (representando las nuevas exposiciones a las dos nuevas cepas) y el espacio de producto \(P * Q\) es el par de lanzamientos independientes.
P = BoxModel([0, 1])
Q = BoxModel([0, 1, 2, 3])
(P * Q).sim(3) Cuando se trata de múltiples variables aleatorias, es común especificar la distribución marginal de cada una y asumir la independencia. Por ejemplo, con \(X\) y \(Y\) con \(X∼Poisson.-1\) y \(Y∼Poisson.-2\), para las variables aleatorias independientes, la distribución conjunta es el producto de las distribuciones marginales.
X, Y = RV(Poisson(1) * Poisson(2))3.3 Condicionamiento
Cuando un resultado \(ω\) se extrae del espacio muestral, un evento ocurre (\(ω∈A\)) o no (\(ω∈A^c\)), así, a simulación de un evento \(A\) devuelve Verdadero para los resultados cuando ocurre el evento y Falso en caso contrario.
Los eventos también se pueden definir y simular. Por razones de programación, los eventos se incluyen entre paréntesis en () lugar de llaves {}.
A = (Y < 3) # evento
A.sim(10)Por ejemplo, para la igualdad lógica, se considera el evento \(\{Y=3\}\)¸.
(Y == 3).sim(10000).tabulate()Como segundo ejemplo, para \(Z=X+Y\) donde \(X\) y \(Y\) son independientes \(X∼Poisson.-1\) y \(Y∼Poisson.-2\), se simula el evento \(A={Z>4}\)
X , Y = RV(Poisson(1) * Poisson(2))
Z = X + Y
A = (Z > 4)
A.sim(10000).tabulate() El condicionamiento de un evento se logra con la barra vertical \(|\). Continuando, se puede simular la distribución condicional de \(X\) dada \({Z = 5}\).
(X | (Z == 5)).sim(10000).plot(normalize = False)3.4 Variables aleatorias múltiples
Para simular variables aleatorias múltiples (parejas) se usa &. También se trabaja con resúmenes y resúmenes normalizados.
xy = (X & Y).sim(100)
xy
xy.tabulate()
xy.tabulate(normalize = True)Las parejas se pueden representar de manera bidimensional en un diagrama de dispersión y se pueden alterar ligeramente para que los puntos no coincidan.
xy.plot(jitter = True)
plt.show()Para dos variables discretas, la gráfica tipo 'tile' produce un mapa de calor donde rectángulos representan las parejas simuladas con sus frecuencias relativas visualizadas en una escala de color.
xy.plot('tile')
plt.show()3.5 Simulación de resultados y variables aleatorias
Llamar (X & Y).sim(10) produce la simulación de los resultados, pero solo (X & Y) se guardan y muestran los valores de. Los resultados de las tiradas se generan en segundo plano pero no se guardan.
Se puede crear una RV que devuelva los resultados del espacio de probabilidad. La función de mapeo predeterminada para RV es la función identidad, \(g(u)=u\), por lo que simular los valores de U = RV(P) devuelve los resultados del BoxModel que P representa. Para simular y mostrar los resultados junto con los valores de \(X\) y \(Y\) usando &.
U = RV(P)
U.sim(10)
(U & X & Y).sim(10)Debido a que el espacio de probabilidad P devuelve pares de valores, U = RV(P) define un vector aleatorio y los componentes individuales de U teambién se pueden desempacar.
U1, U2 = RV(P)
(U1 & U2 & X & Y).sim(10)4 Procesos estocásticos
A continuación se define la idea general de procesos estocásticos y se discutien algunas de sus propiedades a la manera de “Probability!” (2017).
4.1 Definiciones
Proceso estocástico: Una variable aleatoria que cambia con el tiempo.
Se denotan variables aleatorias con letras mayúsculas como \(X\) y variables aleatorias estocásticas con \(X_t\), que simplemente significa el valor de \(X\) en el tiempo \(t\).
El proceso Bernoulli para \(t=1,2,...\) (esto se llama índice del proceso) con todas las variables aleatorias \(X_t\) independientes es un proceso estocástico discreto porque el índice es discreto. Los procesos estocásticos con índices continuos (intervalos continuos de la línea real) se consideran procesos estocásticos continuos.
Este proceso estocástico es solo una serie de ceros y unos que toma el valor \(1\) con probabilidad \(p\) y el valor \(0\) con probabilidad \(1-p\)
# generar proceso Bernoulli
data <- tibble(t = 1:100,
X = rbinom(100, 1, 1 / 2))
# visualizar
ggplot(data, aes(x = t, y = X)) +
geom_line() +
theme_bw() +
ggtitle("Proceso Bernoulli") +
ylim(c(-2, 2))
4.2 Propiedades
Un proceso estocástico es fuertemente estacionario si sus propiedades aleatorias no cambian con el tiempo, o más específicamente si la distribución conjunta de variables aleatorias en diferentes puntos es invariante en el tiempo, i.e.,
\[\begin{equation} \begin{gathered} F_{t_1,t_2,...,t_n}=F_{t_1+k,t_2+k,...,t_n+k}, \end{gathered} \tag{11} \end{equation}\]
Para toda \(k\), \(n\) y conjuntos de intervalos de tiempo \(t_1,t_2,...\) en el proceso estocástico.
Primero, \(F_{t_1}\) es la CDF (función de distribución acumulativa) de \(X_{t_1}\), o el proceso estocástico en el momento \(t_1\) (\(F_{t}=P(X_t<X)\)). Luego, \(F_{t_1,t_2}=P(X_{t_1}<x_1,X_{t_2}<X2)\). Lo que esta ecuación está diciendo, entonces, es que un proceso estocástico es fuertemente estacionario si la distribución conjunta de \(X_t\) a veces \(t_1,t_2,...,t_n\) es igual a la distribución conjunta de \(X_t\) a veces \(t_{1+k},t_{2+k},...,t_{n+k}\) (es decir, los índices de tiempo \(t_1,t_2\) etc. cambiado por \(k\). Por ejemplo, ¿es la distribución conjunta del proceso estocástico en los puntos de tiempo \(1\) y \(2\) la misma que la distribución conjunta del proceso estocástico en los puntos de tiempo \(3\) y \(4\)?.
La estacionariedad débil, naturalmente, tiene condiciones menos restrictivas que la estacionaria fuerte. Un proceso estocástico \(X_t\) es débilmente estacionario si cumple estas tres condiciones:
La media del proceso es constante. Es decir, \(E(X_t)={\mu}\) (donde \({\mu}\) es una constante) para todos los valores de \(t\).
El segundo momento de \(X_t\), o \(X_t^2\), es finito.
Se cumple: \(Cov(X_{t_1},X_{t_2})=Cov(X_{t_{1-k}}, X_{t_{2-k}}).\)
Para todos los conjuntos de \(t_1\) y \(t_2\) y todos los valores de \(k\). Es decir, la covarianza entre dos puntos solo cambia a medida que cambia la distancia entre ellos (es decir, \(k\)), no a medida que avanza en el tiempo . Es decir \(Cov(X_3,X_5)\), debe ser igual a \(Cov(X_{11},X_{14})\) porque \(X_5\) está tan lejos de \(X_8\) como \(X_{14}\) de \(X_{11}\) (tres puntos de tiempo de distancia).
Coloquialmente, un proceso estocástico es una martingala si se puede esperar que la siguiente etapa sea igual a la anterior.
Matemáticamente, un proceso \(X_t\) es martingala si:
\[\begin{equation} \begin{gathered} E(X_{n+1}|X_1,...X_n)=E(X_n). \end{gathered} \tag{12} \end{equation}\]
Para cualquier tiempo \(n\) y con \(E(|X_n|)\) finito.
Así, a cualquier momento \(n\), un proceso estocástico es una martingala si, condicionando en toda la historia observada del proceso hasta tiempo \(n\) (comunmente llamada filtración), la esperanza para la siguiente etapa \(X_{n+1}\) es sólo \(X_n\) (la cual ya es observada, al tiempo \(n\)).
Dos procesos estocásticos \(X_t\) y \(Y_t\) son independientes si, para cualquier secuencia \((t_k,t_{k+1},...,t_{k+n})\) donde \(k\) y \(n\) son números enteros positivos menores que \(\infty\), se tiene que el vector \((X_{t_k},X_{t_{k+1}},...,X_{t_{k+n}})\) es independiente del vector \((Y_{t_k},Y_{t_{k+1}},...,Y_{t_{k+n}})\).
Coloquialmente, para que \(X_t\) y \(Y_t\) sean independientes, cualquier subconjunto de \(X_t\) debe ser independiente del mismo subconjunto (en los mismos puntos de tiempo) de \(Y_t\).
\(X_t\) tiene densidad \(f_{x_t}=(x_{t_{k}},x_{t_{k+1}},...,x_{t_{k+n}})\), donde \(f_{x_t}\) se define simplemente como la densidad conjunta de \(X_t\) (si cada \(X_t\) fuera i.i.d., \(f_{x_t}\) es sólo el producto de las densidades marginales de \(X_t\) en todos los puntos de tiempo individuales). Similarmente, la densidad conjunta de \(Y_t\) tiene densidad \(f_{y_t}=(y_{t_{k}},y_{t_{k+1}},...,y_{t_{k+n}})\).
Los vectores son independientes si la densidad conjunta de ambos vectores \(f_{x_t,y_t}\) es igual al producto de \(f_{x_t}\) y \(f_{y_t}\): \(f_{x_t,y_t}(x_{t_{k}},x_{t_{k+1}},...,x_{t_{k+n}},y_{t_{k}},y_{t_{k+1}},...,y_{t_{k+n}})=f_{x_t}(x_{t_{k}},x_{t_{k+1}},...,x_{t_{k+n}})\times \\\) \(f_{y_t}(y_{t_{k}},y_{t_{k+1}},...,y_{t_{k+n}})\).
Coloquialmente, un proceso estocástico con la propiedad de Markov solo se preocupa por su estado más reciente. Esto es, el proceso estocástico \(X_t\) satisface la propiedad de Markov si:
\[\begin{equation} \begin{gathered} P(X_t|X_1,X_2,...,X_{t−1})=P(X_t|X_{t−1}). \end{gathered} \tag{13} \end{equation}\]
Es decir, dentro de la filtración de \(X_t\) (toda la información antes del tiempo \(t\)), dice toda la información sobre \(X_t\) (\(X_t\) es condicionalmente independiente de \(X_1,X_2,...,X_{t-2}\) dado \(X_{t-1}\)).
Sea \(X_t\) un proceso estocástico continuo tal que, para todos \(t\) y \(k<t\), se tiene que \(X_t-X_{t-k}\) es una variable aleatoria simétrica alrededor de \(0\) (una v.a. es simétrica alrededor del \(0\) si \(Y\) tiene la misma distribución que \(−Y\)). Esto es, el incremento del proceso (cantidad agregada al proceso en cada intervalo) es simétrico.
5 Aplicaciones
En los modelos estocásticos, al ser modelos en los que las variables cambian aleatoriamente, es necesario ejecutar el modelo no solo una sino muchas veces y luego estudiar la distribución resultante de los resultados. Al mismo tiempo, la mayoría de los procesos biológicos son intrínsecamente estocásticos, por lo que los modelos estocásticos suelen ser más realistas que los deterministas.
Las apĺicaciones se encuentran en ecología (por ejemplo, modelos neutrales de biodiversidad), genética de poblaciones (por ejemplo, el proceso de coalescencia) y evolución (por ejemplo, métodos para estimar árboles filogenéticos). Una amplia comprensión de clases importantes de modelos estocásticos, nuevamente con ejemplos de la ecología y la evolución, puede resumirse en la Tabla 5 (Engelstädter 2019).
| Proceso estocástico | Tiempo | Caracteristicas | Ejemplos de | |
|---|---|---|---|---|
| Proceso de ramificación | discreto o continuo | modelo de población donde el número de descendientes de cada individuo se extrae de la misma distribución | colonización de nuevos hábitats, propagación de nuevas enfermedades | |
| Cadena de Markov | discreto | cambia entre diferentes estados, con probabilidades que dependen del estado anterior | Sustituciones de nucleótidos en la evolución de la secuencia de ADN. | |
| Proceso Poisson | continuo | eventos que suceden de forma independiente y con una pequeña probabilidad por unidad de tiempo | mutaciones en el linaje de los individuos, proceso coalescente | |
| Movimiento Browniano | continuo | cambios aleatorios en la variable pero con media cero | movimiento de individuos en el espacio |
Todos estos procesos de modelos pueden considerarse casos especiales de una clase más amplia de procesos estocásticos denominados procesos de Markov. Los procesos de Markov se caracterizan por ser sin memoria, lo que significa que el cambio en las variables está influenciado solo por su presente pero no por sus estados pasados.
5.1 Procesos de ramificación
Los procesos de ramificación pueden entenderse como modelos de población en los que el número de descendientes son variables aleatorias independientes y distribuidas de forma idéntica. Por ejemplo, supongamos una población en la que cada individuo \(i\) produce descendencia \(k_i\), con \(k_i\) de \(0\) a \(2\). La distribución de \(k_i\) está dada de acuerdo con la siguiente función:
\[\begin{equation} \begin{gathered} k_i=\begin{cases} 0 & \text{con probabilidad} \,\ p_0=0.2 \\ 1 & \text{con probabilidad} \,\ p_1=0.5 \\ 2 & \text{con probabilidad} \,\ p_2=0.3 \end{cases} \end{gathered} \tag{14} \end{equation}\]
Denotando el tamaño de la población en el tiempo \(t\) por \(N_t\) y asumiendo que todos los individuos adultos mueren después de la reproducción, se describe al modelo como
\[\begin{equation} \begin{gathered} N_t=\begin{cases} 0 & \text{si} \,\ N_t=0 \\ \sum_{i=1}^N k_i & \text{si} \,\ N_t>0 \end{cases} \end{gathered} \tag{15} \end{equation}\]
Dada la dinámica poblacional de este modelo, no es difícil calcular las probabilidades de la ecuación (14) que el número esperado de descendientes es \(1.1\) para todos los individuos (\(E[k_i]=0\times 0.2+1\times 0.5+2\times 0.3=1.1\)). Por lo tanto, en promedio, se espera que la población pudiera crecer. Sin embargo, como se trata de un modelo estocástico, la extinción es posible y de hecho, una cuestión importante en el proceso de ramificación es con qué probabilidad se produce la extinción. Ésta y otras preguntas se pueden abordar utilizando funciones generadoras de probabilidad.
5.2 Cadenas de Markov
Los procesos de ramificación son un caso especial de una clase de procesos estocásticos de tiempo discreto llamados cadenas de Markov. En una cadena de Markov, una variable puede tomar varios valores discretos que se denominan estados. En cada paso de tiempo, esos estados pueden cambiar a otros estados con ciertas probabilidades.
Como ejemplo, sea la dinámica de las sustituciones de nucleótidos en un sitio determinado en el ADN de una especie. Los cuatro estados son los cuatro nucleótidos, es decir, la variable de nuestro modelo puede tomar cualquiera de los cuatro valores A, C, T y G. Mediante mutación y deriva, este nucleótido puede ser reemplazado por un nucleótido diferente dentro de un cierto período de tiempo. En principio, se podría asignar cada una de las doce sustituciones posibles (A→ C, A→ G, A →T, C →A, C →G, etc.) su propia probabilidad.
Las probabilidades de transición entre diferentes purinas o pirimidinas (A <-> G o C <-> T) ocurren con probabilidad \(α\) mientras que las transversiones (de purina a pirimidina o viceversa) ocurren con una probabilidad menor \(β\). Por lo tanto, en total, la probabilidad de que una base cambie en un paso de tiempo es \(α+2β\) y la probabilidad de que no cambie es \(1-α-2β\).
No hay un estado absorbente en este modelo, es decir, ningún estado en el que el sistema se atasque y ya no cambie. Además, el modelo es completamente simétrico: aunque al principio permanece en el nucleótido original durante algún tiempo, a la larga, los cuatro estados tienen las mismas posibilidades de alcanzarse. En otras palabras, la distribución de probabilidad de equilibrio de los cuatro estados es el vector ( \(\frac{1}{4}\),\(\frac{1}{4}\),\(\frac{1}{4}\),\(\frac{1}{4}\)) en este modelo.
Como segundo ejemplo, se considera un modelo de extinción de especies en el que una especie puede tener tres estados de conservación diferentes: sin preocupación (N), amenazada (T) y extinta (E). Suponiendo que dentro de un paso de tiempo una especie solo puede cambiar de una de estas categorías a una adyacente, las especies pueden pasar de N a T a una velocidad determinada por α, y de T a N a una velocidad determinada por β. Una vez que las especies pasan de T a E (a una velocidad γ), permanecerán en este estado de forma permanente. El estado absorbente es el estado E. En este modelo, cuando α>0 y γ>0, la etapa E eventualmente se alcanzará y luego se mantendrá indefinidamente.
Las cadenas de Markov y sus diversas extensiones y aplicaciones juegan un papel central en la biología computacional. Una extensión importante son los modelos ocultos de Markov. Aquí, los estados del proceso están ocultos, es decir, no se pueden observar. En cambio, hay otros estados observables que son generados ( emitidos ) con ciertas probabilidades por los estados ocultos. Los modelos ocultos de Markov se utilizan ampliamente en bioinformática, por ejemplo, en métodos de predicción de plegamiento de proteínas, alineación de secuencias de ADN y predicción de genes.
Una aplicación importante de las cadenas de Markov son los métodos de Markov Chain Monte Carlo (MCMC). Hay muchos problemas en biología en que hay que estimar distribuciones de probabilidad complejas. Un ejemplo que surge en el campo de la filogenética es la distribución de probabilidades de diferentes árboles filogenéticos que explican la relación evolutiva entre especies, dados ciertos datos de secuencia de ADN para la especie. Dentro de un marco bayesiano, se puede estimar la probabilidad de un árbol dados los datos (la probabilidad posterior) a partir de la probabilidad de observar estos datos dado el árbol (la probabilidad del árbol). Esto se logra mediante varios modelos de evolución de la secuencia de ADN. Desafortunadamente, debido a que el número de árboles posibles suele ser descomunal, es imposible calcular esta probabilidad para todos los árboles. Aquí es donde entra MCMC: en lugar de intentar atravesar todos los árboles (imposible) o seleccionar un conjunto aleatorio de árboles (muy ineficiente), se realiza una especie de paseo aleatorio por el bosque de árboles. Comenzando con un árbol elegido al azar, se estima su probabilidad. Luego, se elige aleatoriamente otro árbol que es similar al primero y también se determina su probabilidad. Dependiendo de estas dos probabilidades y del algoritmo exacto que se utilice (una opción popular es el algoritmo Metropolis-Hastings), se procede al nuevo árbol o se queda con el antiguo. Este proceso se repite muchas veces. El número relativo de pasos de tiempo que el proceso pasa en los diferentes estados en MCMC convergerá a la distribución de probabilidad subyacente de estos estados. En otras palabras, MCMC permite obtener muestras de una distribución de probabilidad compleja.
5.3 Proceso Poisson
Los procesos de Poisson son procesos en los que ciertos eventos ocurren independientemente entre sí y con una probabilidad pequeña y constante por unidad de tiempo. Por ejemplo, estos eventos podrían ser clientes que llegan a una tienda, partículas radiactivas en descomposición o estrellas fugaces. El procesos Poisson se ha estudiado durante mucho tiempo y tiene una serie de propiedades simples que lo hace atractivo. Más importante aún, el número de eventos que ocurren dentro de un intervalo de tiempo dado sigue una distribución Poisson. Además, el intervalo de tiempo entre dos eventos sucesivos sigue una distribución exponencial.
En ecología teórica y biología evolutiva, los procesos Poisson tienen un papel central. Por ejemplo, al considerar una especie en evolución a lo largo del tiempo, la acumulación de mutaciones en esta especie puede modelarse como un proceso Poisson y ser la base para estimar la velocidad a la que se acumulan las mutaciones, o emplearse para estimar el momento en el que una especie ancestral se dividió en dos especies existentes. Los procesos Poisson también pueden ser parte de procesos estocásticos más complejos como los procesos de coalescencia (para aplicaciones modernas en ecología de poblaciones).
5.4 Movimiento Browniano
Son procesos de tiempo continuo en los que una o varias variables están sujetas a pequeñas fluctuaciones aleatorias. Más precisamente, los cambios en la (s) variable (s) en pasos de tiempo infinitesimalmente pequeños siguen una distribución normal con media cero y una varianza constante que determina la velocidad del proceso. Cualquier cambio en la (s) variable (s) es independiente de los cambios previos (la propiedad de Markov o falta de memoria). El resultado es una especie de caminata aleatoria donde tanto el tiempo como el espacio son continuos.
Se pueden utilizar, por ejemplo, como modelos nulos para el movimiento de organismos en el espacio. Tal movimiento podría describirse como una especie de deambular sin rumbo fijo, en contraste quizás con una exploración más sistemática del hábitat o un movimiento dirigido hacia algún objetivo. El Browniano también se ha utilizado ampliamente para modelar cómo cambian los rasgos continuos a lo largo del tiempo en un clado de especies en evolución.
Session info
## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
##
## Matrix products: default
## BLAS: /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3
## LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/liblapack.so.3
##
## locale:
## [1] LC_CTYPE=es_MX.UTF-8 LC_NUMERIC=C
## [3] LC_TIME=es_MX.UTF-8 LC_COLLATE=es_MX.UTF-8
## [5] LC_MONETARY=es_MX.UTF-8 LC_MESSAGES=es_MX.UTF-8
## [7] LC_PAPER=es_MX.UTF-8 LC_NAME=C
## [9] LC_ADDRESS=C LC_TELEPHONE=C
## [11] LC_MEASUREMENT=es_MX.UTF-8 LC_IDENTIFICATION=C
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] RColorBrewer_1.1-2 kableExtra_1.3.4 knitr_1.33 EpiCompare_0.0.1
## [5] ggtern_3.3.5 forcats_0.5.1 stringr_1.4.0 dplyr_1.0.7
## [9] purrr_0.3.4 readr_2.0.0 tidyr_1.1.3 tibble_3.1.3
## [13] ggplot2_3.3.5 tidyverse_1.3.1 BiocStyle_2.20.2
##
## loaded via a namespace (and not attached):
## [1] fs_1.5.0 lubridate_1.7.10 webshot_0.5.2
## [4] httr_1.4.2 latex2exp_0.5.0 tensorA_0.36.2
## [7] tools_4.1.1 backports_1.2.1 bslib_0.2.5.1
## [10] utf8_1.2.2 R6_2.5.0 DBI_1.1.1
## [13] colorspace_2.0-2 withr_2.4.2 tidyselect_1.1.1
## [16] gridExtra_2.3 bayesm_3.1-4 compiler_4.1.1
## [19] compositions_2.0-2 cli_3.0.1 rvest_1.0.0
## [22] xml2_1.3.2 labeling_0.4.2 bookdown_0.22
## [25] sass_0.4.0 scales_1.1.1 DEoptimR_1.0-9
## [28] robustbase_0.93-9 systemfonts_1.0.2 digest_0.6.27
## [31] rmarkdown_2.9 svglite_2.0.0 pkgconfig_2.0.3
## [34] htmltools_0.5.2 dbplyr_2.1.1 fastmap_1.1.0
## [37] highr_0.9 rlang_0.4.11 readxl_1.3.1
## [40] rstudioapi_0.13 jquerylib_0.1.4 generics_0.1.0
## [43] farver_2.1.0 jsonlite_1.7.2 magrittr_2.0.1
## [46] Matrix_1.3-4 Rcpp_1.0.7 munsell_0.5.0
## [49] fansi_0.5.0 reticulate_1.20 viridis_0.6.1
## [52] proto_1.0.0 lifecycle_1.0.1 stringi_1.7.4
## [55] yaml_2.2.1 MASS_7.3-54 plyr_1.8.6
## [58] grid_4.1.1 crayon_1.4.1 lattice_0.20-45
## [61] haven_2.4.1 hms_1.1.0 pillar_1.6.1
## [64] reprex_2.0.0 glue_1.4.2 evaluate_0.14
## [67] BiocManager_1.30.16 modelr_0.1.8 png_0.1-7
## [70] vctrs_0.3.8 tzdb_0.1.2 cellranger_1.1.0
## [73] gtable_0.3.0 assertthat_0.2.1 xfun_0.26
## [76] broom_0.7.8 viridisLite_0.4.0 ellipsis_0.3.2