Licenciado bajo Creative Commons BY-NC-SA 4.0.
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
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)
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)