Licenciado bajo Creative Commons BY-NC-SA 4.0.

1 Cadenas de Markov de tiempo discreto

1.1 Introducción

Para la manipulación de cadenas de Markov en R, vamos a utilizar la librería markovchain, creada por Giorgio Alfredo Spedicato. A continuación, nos familiarizaremos con su uso a través del siguiente modelo simple de predicción del tiempo visto en teoría (C: calor, F: frío):

\[\xymatrix{ *=<5mm,5mm>[o][F]{C} \ar@(ul,dl)[]_{0.6} \ar@/^/[r]^{0.4} & *=<5mm,5mm>[o][F]{F} \ar@(dr,ur)[]_{0.7} \ar@/^/[l]^{0.3} }\qquad\qquad P = \begin{pmatrix} 0.6 & 0.4 \\ 0.3 & 0.7 \end{pmatrix}\]

Esta introducción está basada en la documentación oficial.

En primer lugar, debemos importar la librería para tener acceso a sus clases y métodos. La siguiente instrucción equivale al #include de C o el import de Python:

library(markovchain)
## Warning: package 'markovchain' was built under R version 3.6.2

1.2 Matriz de transición

Definimos la matriz de transición del ejemplo de arriba. Una matriz en R se define mediante el método matrix (véase la ayuda). Como entrada, recibe un vector de datos unidimensional, por lo que hay que especificar al menos una dimensión (número de filas o columnas) y cómo debe interpretar dicho vector (por filas o por columnas). En este caso, introducimos una matriz de 2x2 (ncol = 2) por filas (byrow = TRUE):

P <- matrix(data = c(0.6, 0.4,
                     0.3, 0.7), 
            ncol = 2, byrow = TRUE)

1.3 Cadenas de Markov en R

Una cadena de Markov vendrá representada por un objeto de tipo “markovchain”. Para su creación, utilizamos el método genérico new e introducimos los nombres de los estados y la matriz de transición. Opcionalmente, podemos darle un nombre a la cadena, como se ve a continuación.

weather <- new("markovchain", name = "Simple weather model", 
               states = c("C", "F"), transitionMatrix = P)
weather
## Simple weather model 
##  A  2 - dimensional discrete Markov Chain defined by the following states: 
##  C, F 
##  The transition matrix  (by rows)  is defined as follows: 
##     C   F
## C 0.6 0.4
## F 0.3 0.7
summary(weather)
## Simple weather model  Markov chain that is composed by: 
## Closed classes: 
## C F 
## Recurrent classes: 
## {C,F}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE
plot(weather)

En R no podemos utilizar los operadores habituales de multiplicación * y exponenciación ^ con matrices. La multiplicación se realiza con el siguiente operador:

P %*% P
##      [,1] [,2]
## [1,] 0.48 0.52
## [2,] 0.39 0.61

Y para elevar una matriz a una potencia, necesitamos cargar la librería expm:

library(expm)
## Loading required package: Matrix
## 
## Attaching package: 'expm'
## The following object is masked from 'package:Matrix':
## 
##     expm
P%^%4
##        [,1]   [,2]
## [1,] 0.4332 0.5668
## [2,] 0.4251 0.5749

En cambio, el objeto markovchain se comporta aritméticamente como una matriz de transición sobre la que sí podemos aplicar los operadores habituales. La siguiente tabla detalla las principales operaciones definidas sobre él.

Método Resultado
mcA * mcB Multiplicación de matrices de transición
mcA^4 Exponenciación de la matriz de transición
mcA[1,] Acceso a la fila 1 de la matriz de transición
mcA[, 3] Acceso a la columna 3 de la matriz de transición
mcA[3, 2] Acceso al elemento (3, 2) de la matriz de transición
t(mcA) Trasposición de la matriz de transición
mcA == mcB Comparación entre dos objetos (estados y matriz)
dim(mcA) Dimensión de la matriz de transición
states(mcA) Extracción de los estados

Por tanto, podemos aplicar este objeto directamente a la ecuación de Chapman-Kolmogorov:

\[\pi^{(n)} = \pi^{(0)}P^n\]

Partimos del estado inicial “F” y obtenemos el vector de distribución de probabilidad para \(n = 6\):

initialState <- c(0, 1)
newState <- initialState * weather^6
newState
##             C        F
## [1,] 0.428259 0.571741

Pero, además, el objeto markovchain nos permite obtener mucha información sobre la cadena: irreducibilidad, accesibilidad, clases, tipos de estados, etc.

is.irreducible(weather)
## [1] TRUE
is.accessible(weather, from = "C", to = "F")
## [1] TRUE
period(weather)
## [1] 1
absorbingStates(weather)
## character(0)
transientStates(weather)
## character(0)
steadyStates(weather)
##              C         F
## [1,] 0.4285714 0.5714286

Lo más interesante de esta librería es que nos permite realizar de forma sencilla simulación y estimación de cadenas de Markov. La simulación, dada una cadena de Markov previamente definida, se realiza de manera análoga a la simulación de variables aleatorias convencionales (rnorm, rexp, etc.) mediante la función rmarkovchain.

Generación de 1000 muestras aleatorias a partir de la cadena de Markov que hemos definido al comienzo del tutorial. El estado inicial es aleatorio:

x <- rmarkovchain(1000, weather)

Generación de 1000 muestras aleatorias partiendo del estado inicial “C”:

y <- rmarkovchain(1000, weather, t0 = "C")

La estimación consiste en ajustar unas probabilidades de transición a partir de un vector de muestras aleatorias. El ajuste a partir del vector x generado anteriormente se aproxima bastante a la cadena original:

weatherFit <- markovchainFit(x)
plot(weatherFit$estimate)

1.4 Convergencia

Dada una cadena de Markov irreducible y aperiódica, se puede demostrar que existe un único estado estacionario al que la cadena converge cuando \(n\to\infty\) independientemente del estado inicial. Tomando nuestro modelo de predicción del tiempo, podemos ilustrar esta propiedad elevando la matriz de transición a una potencia suficientemente alta y observando que todas las filas de la matriz resultante tienden al estado estacionario:

is.irreducible(weather)
## [1] TRUE
weather^100
## Simple weather model^100 
##  A  2 - dimensional discrete Markov Chain defined by the following states: 
##  C, F 
##  The transition matrix  (by rows)  is defined as follows: 
##           C         F
## C 0.4285714 0.5714286
## F 0.4285714 0.5714286
steadyStates(weather)
##              C         F
## [1,] 0.4285714 0.5714286

Esto es equivalente a decir que, cuando simulamos una cadena de Markov (obtenemos un conjunto de muestras), las probabilidades de los estados serán estables al alcanzar un número de muestras lo suficientemente elevado, y no van a cambiar por añadir más.

samples <- rmarkovchain(1000, weather)

# Proporción de cada estado acumulada en cada paso
yC <- cumsum(samples=="C") / seq_along(samples)
yF <- cumsum(samples=="F") / seq_along(samples)
# Evolución de la distribución de C
plot(yC, type="l", col="red", ylim=c(0,1))
# Evolución de la distribución de F
lines(yF, col="green")
# Valores estacionarios teóricos
abline(steadyStates(weather)[1], 0, lty=2, col="red")
abline(steadyStates(weather)[2], 0, lty=2, col="green")

En la figura anterior se aprecia que la simulación converge. Pero, ¿qué potencia es lo suficientemente alta, qué número de muestras el lo suficientemente elevado para asegurar una buena aproximación al estado estacionario? Aquí entra en juego el llamado tiempo de convergencia, que nos dará una idea de a qué velocidad se va a reducir el error conforme avanza el tiempo.

Podemos definir el error, de forma simple, como el sumatorio de las diferencias absolutas entre componentes de dos potencias sucesivas de la matriz de transición, es decir,

\[\epsilon(n) = \sum_i\sum_j|\left(P^n\right)_{ij} - \left(P^{n-1}\right)_{ij}|\]

Vamos a construir una función en R para obtener este error y poder visualizarlo.

# Función de error
# @param n es el tiempo (valor o conjunto de valores)
# @param mc es la cadena de Markov
err <- function(n, mc) {
  # Comprobamos que los valores de entrada son números naturales
  if (!isTRUE(all(n > 0 && n == floor(n)))) 
    stop("'n' must only contain positive integer values")
  
  # Reservamos memoria para la salida
  res <- numeric(length(n))
  # Para cada valor de entrada "n", calculamos err(n)
  for (i in 1:length(n)) {
    Pn <- (mc^n[i])@transitionMatrix
    Pn1 <- (mc^(n[i]-1))@transitionMatrix
    res[i] <- sum(abs(Pn - Pn1))
  }
  return(res)
}

# Representación de la función de error
x <- 1:10
y <- err(x, weather)
plot(x, y, type="o")

Comprobamos que el error decae exponencialmente.

2 Cadenas de Markov de tiempo continuo

Ejemplo 1: Sea una gasolinera con un solo surtidor y sin espacio para que los vehículos esperen (si un vehículo llega a la gasolinera y el surtidor está ocupado, se va). A la gasolinera llegan vehículos según un proceso de Poisson de tasa \(\lambda\) igual a 3 cada 20 minutos, de los cuales \(prob(c)=\) 75% son coches y \(prob(m)=\) 25% motocicletas. El tiempo que le lleva a un usuario dejar libre el surtidor (es decir, repostar y pagar) se puede modelar como una variable aleatoria exponencial, de media \(1/\mu_c = 8\) minutos si es un coche, y \(1/\mu_m=3\) minutos si es una moto.

Este problema se puede modelar con la siguiente cadena de Markov de tiempo continuo:

\[\xymatrix{ *=<15mm,8mm>[o][F]{car} \ar@/^/[r]^{\mu_c} & *=<15mm,8mm>[o][F]{empty} \ar@/^/[l]^{\lambda \cdot prob(c)} \ar@/^/[r]^{\lambda \cdot prob(m)} & *=<15mm,8mm>[o][F]{m/cycle} \ar@/^/[l]^{\mu_m} }\qquad\qquad Q = \begin{pmatrix} -1/8 & 1/8 & 0 \\ 0.75\cdot 3/20 & -3/20 & 0.25\cdot 3/20 \\ 0 & 1/3 & -1/3 \end{pmatrix}\]

La cadena es irreducible y recurrente. Para encontrar la distribución de estado estacionaria teóricamente, debemos resolver las ecuaciones de balance

\[pQ = 0\]

a las que hay que añadir la relación de cierre

\[\sum_i p_i = 1\]

para obtener un sistema compatible determinado. Dado que hay \(\operatorname{dim}(Q)-1\) columnas independientes, dicha relación de cierre se traduce en sustituir cualquier columna por unos e igualar a uno dicha componente al otro lado de la ecuación, es decir:

\[p\begin{pmatrix} 1 & 1/8 & 0 \\ 1 & -3/20 & 0.25\cdot 3/20 \\ 1 & 1/3 & -1/3 \end{pmatrix} = (1, 0, 0)\]

Al resolver el sistema, obtendremos un vector \(p\) que representa la probabilidad de estar en cada uno de los estados a largo plazo. Por tanto, para calcular el número medio de usuarios, estamos interesados en sumar las probabilidades multiplicadas por el número de usuarios en sus respectivos estados. En este caso, \(1\cdot p_1 + 0\cdot p_2 + 1\cdot p_3\).

lambda <- 3/20    # Tasa de llegadas
mu <- c(1/8, 1/3) # Tasa de servicio de coches y motos respectivamente
p <- 0.75         # Probabilidad de que el vehículo sea coche

A <- matrix(c(1,   mu[1],            0,
              1, -lambda, (1-p)*lambda,
              1,   mu[2],       -mu[2]), byrow=T, ncol=3)

P <- solve(t(A), c(1, 0, 0))
N_average_theor <- sum(P * c(1, 0, 1))
N_average_theor
## [1] 0.5031056

En el código anterior, la función t traspone la matriz para resolver el sistema con solve. Consulte la ayuda para más detalles.

3 Resolución de problemas con R

3.1 Problema 5.10

El problema del gato y el ratón. Hay cinco cajas en fila, con un gato en la primera y un ratón en la última. En cada intervalo de tiempo \(T\), cada animal salta a una caja colindante, completamente al azar. El juego termina cuando ambos animales saltan a la misma caja. Calcule la duración media del juego.

3.2 Problema 5.5

Una tarjeta de red tiene tres velocidades de transmisión, 1, 4 y 8 Mbps, que dan lugar a una probabilidad de pérdida de trama de \(p_1=1/2\), \(p_4=1/2\) y \(p_8=1/4\) respectivamente. El esquema que se usa para adaptar la velocidad de transmisión es el siguiente: ante una pérdida, se transmite a 1 Mbps, mientras que tras dos éxitos consecutivos a una velocidad se aumenta la velocidad al siguiente valor superior disponible. Calcule la velocidad media (en Mbps) de transmisión.

3.3 Problema 6.8

Sea un edificio con dos ascensores, cada uno con un tiempo de vida que se puede modelar con una variable aleatoria exponencial de media \(1/\lambda\). Hay dos políticas de reparación:

  • Se avisa a mantenimiento cuando los dos ascensores están estropeados, siendo el tiempo de reparación de los dos ascensores (a la vez) una variable aleatoria exponencial de media \(1/\mu\).
  • Se avisa a mantenimiento en cuanto un ascensor se estropea, siendo el tiempo de reparación por ascensor una variable aleatoria exponencial de media \(2/\mu\), y sólo se repara un ascensor a cada vez.

Calcule la proporción de tiempo que los dos ascensores están funcionando, si \(1/\lambda = 2\) semanas y \(1/\mu= 1\) semana.