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)