Cadenas de Markov en R

Procesos Estocásticos

Instalación del paquete

Para la instalación de este paquete en R se utiliza el comando “install.package()” y para el estudio de cadenas de Markov vamos a utilizar el paquete “markovchain” en su version 0.6.9.11 publicado por : Giorgio Alfredo Spedicato, Tae Seung Kang y otros.

# Comando para instalar el paquete
# install.packages("markovchain")
library(markovchain)

## Package:  markovchain
## Version:  0.8.5-4
## Date:     2021-01-07
## BugReport: https://github.com/spedygiorgio/markovchain/issues

Introducción a markovchain

Esta libreria pretende proveer objetos para realizar análisis estadísticos de Cadenas de markov a tiempos discretos. Asumamos que tenemos una cadena de markov \(X=\{X_1 ~,X_2,...\}\) definida en el espacio de estados \(S = \{a,b,c\}\) cuya matriz de transición es: \[P=(\left( {\begin{array}{ccc} 0 & 0.5 & 0.5 \\ 0.5 & 0 & 0.5 \\ 0.5 & 0.5 & 0 \\ \end{array} } \right)\]

Dicha cadema podemos crearla en R, de la siguiente forma:

Crear la matriz de transición P:

P = matrix(c(0,0.5,0.5,0.5,0,0.5,0.5,0.5,0),nrow = 3,byrow = TRUE)
P

##      [,1] [,2] [,3]
## [1,]  0.0  0.5  0.5
## [2,]  0.5  0.0  0.5
## [3,]  0.5  0.5  0.0

El argumento “nrow”es la función a la cual indicamos cuantas filas queremos que tenga nuestra matriz P tenga, y “byrow” es para que almacene los elementoss de la matriz almacenados en c(), fila por fila.

Creamos el objeto markovchain

mc = new("markovchain",transitionMatrix=P,states=c("a","b","c"),name="Cadena 1")

Para hacer una revisión previa al análisis se puede hacer con el comando “str()” y “summary” que devuelve la estructura del objeto

str(mc)

## Formal class 'markovchain' [package "markovchain"] with 4 slots
##   ..@ states          : chr [1:3] "a" "b" "c"
##   ..@ byrow           : logi TRUE
##   ..@ transitionMatrix: num [1:3, 1:3] 0 0.5 0.5 0.5 0 0.5 0.5 0.5 0
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:3] "a" "b" "c"
##   .. .. ..$ : chr [1:3] "a" "b" "c"
##   ..@ name            : chr "Cadena 1"

mc es un objeto de tipo markovchain, que representa una lista con arguentos provisto declarados en el contructor.

El resumen de la cadena es:

summary(mc)

## Cadena 1  Markov chain that is composed by: 
## Closed classes: 
## a b c 
## Recurrent classes: 
## {a,b,c}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE

Para hacer una visualización del diagrama de transición de la cadena lo podemos realizar mediante la función plot

plot(mc)

Otras funciones inportantes son:

absorbingStates(): Identifica los estados Absorbentes
transientStates(): Identifica los estados Transitorios
recurrentClasses(): Identifica las clases Recurrentes

Para la cadena de markov definida se obtiene que:

recurrentStates(mc)

## [1] "a" "b" "c"

transientStates(mc)

## character(0)

recurrentClasses(mc)

## [[1]]
## [1] "a" "b" "c"

Análisis Probabilístico

Para conocer la probabilidad de transición en 1 paso entre un estado y otro basta con utilizar la función transitionProbability(), con los argumentos:

object: la cadena de markov
t0: el estado en el tiempo 0
t1: el estado en el tiempo 1

La probabilida de transición en un paso dfel estado “a” al estado “c” es:

transitionProbability(object=mc, t0="a",t1="c")

## [1] 0.5

Hay que recordar que dicha probabilida es un elemento de la matriz de transición P, por lo tanto, la probabilida de transición del estado “a” al estado “c” es simplemente \(P_{23}\)

mc[2,3]

## [1] 0.5

Es posible computar la matriz de transición en n pasos, simplemente computando la n-ésima potencia de la matriz de transición P, como ejemplo calcularemos la matriz de transición en n=5 pasos.

n=5
mc ^ n

## Cadena 1^5 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  a, b, c 
##  The transition matrix  (by rows)  is defined as follows: 
##         a       b       c
## a 0.31250 0.34375 0.34375
## b 0.34375 0.31250 0.34375
## c 0.34375 0.34375 0.31250

Tambien se puede conocer la distribución de la cadena en n pasos adelante \((P(X_n))\) multiplicando la distribucion inicial de \(X_0\) por la matriz de transición en n pasos \((P^n)\), calcule la distribucion de la cadena en el tiempo n = 6, si la distribucion inicial de la cadena es “(0.5,0.2,0.3)”.

X0=c(0.5,0.2,0.3)
n=6
Xn = X0*(mc^n)

por lo tanto la distribución de la cadena en 6 pasos es:

Xn

##              a       b         c
## [1,] 0.3359375 0.33125 0.3328125

Puesto que \(Xn\) es una función de densidad, la suma de las probabilidades en todos los estados debe ser 1.

sum(Xn)

## [1] 1

Finalmente encontrar la distribucion estacionaria de la cadena se obtiene mediante la función steadyStates de la siguiente forma:

DistEst= steadyStates(mc)
DistEst

##              a         b         c
## [1,] 0.3333333 0.3333333 0.3333333

Hay que recordad que los tiempos medio de recurrencia son los inversos multiplicativos de la distribución estacionaria y puede ser computados facilmente

M = 1/DistEst
M

##      a b c
## [1,] 3 3 3

Ejercicio

Dibuje el diagrama de transición, determine las clases de comunicación de las siguientes cadenas de Markov, clasifique éstas como recurrentes o transitorias, y encuentre la distribución estacionaria si existe.

\[P=\left( {\begin{array}{ccc} 1/2 & 1/2 & 0 & 0 \\ 0 & 1/2 & 1/2 & 0\\ 0 & 1/2 & 1/2 & 0\\ 1/4 & 1/4 & 1/4 & 1/4\\ \end{array} } \right)\]

K = matrix(c(0.5,0.5,0,0,0,0.5,0.5,0,0,0.5,0.5,0,0.25,0.25,0.25,0.25),nrow = 4,byrow = TRUE)
K

##      [,1] [,2] [,3] [,4]
## [1,] 0.50 0.50 0.00 0.00
## [2,] 0.00 0.50 0.50 0.00
## [3,] 0.00 0.50 0.50 0.00
## [4,] 0.25 0.25 0.25 0.25

ejemplo = new("markovchain",transitionMatrix = K,states=c("0","1","2","3"))
summary(ejemplo)

## Unnamed Markov chain  Markov chain that is composed by: 
## Closed classes: 
## 1 2 
## Recurrent classes: 
## {1,2}
## Transient classes: 
## {0},{3}
## The Markov chain is not irreducible 
## The absorbing states are: NONE

plot(ejemplo)

Distribución estacionaria de la matriz

steadyStates(ejemplo)

##      0   1   2 3
## [1,] 0 0.5 0.5 0

Aplicaciones

En aplicaciones cotidianas la matriz de transición P es desconocida, usualmente la información conocida es un registro en el tiempo del comportamiento de algún fenómeno, y por lo tanto es necesario estimar la matriz de transición para los registros dados.

Como ejemplo se presenta el registro semanal del número de veces que llueve en alguna ciudad desconocida, almacenados en la base de datos rain. ## Cargando y revisando los archivos

data(rain)
str(rain)

## 'data.frame':    1096 obs. of  2 variables:
##  $ V1  : int  3 2 2 2 2 2 2 3 3 3 ...
##  $ rain: chr  "6+" "1-5" "1-5" "1-5" ...

La variable rain de los datos guarda el registro del número de veces que llueve, la tabla de frecuencias para rain es:

table(rain$rain)

## 
##   0 1-5  6+ 
## 548 295 253

Podemos observar los primeros registros de los datos de la siguiente forma:

head(rain)

##   V1 rain
## 1  3   6+
## 2  2  1-5
## 3  2  1-5
## 4  2  1-5
## 5  2  1-5
## 6  2  1-5

Como modelo se propone una cadena de markov finita de 3 estados para modelar los registros semanales del numero de veces que llueve almacenados en rain.

PE = rain$rain
head(PE)

## [1] "6+"  "1-5" "1-5" "1-5" "1-5" "1-5"

La libreria markovchain, nos provee herramientas para estimar la matriz de transición de la cadena mediante los datos obtenidos. La función “CreateSequenceMatrix()” crea una matriz de secuencia.

P1 = createSequenceMatrix(PE)
P1

##       0 1-5  6+
## 0   362 126  60
## 1-5 136  90  68
## 6+   50  79 124

Donde cada elemento \(p_{ji}\) de la matriz representa el número de veces que el proceso paso del estado i al estado j. La funcion markovchainFit() estima la matriz de transición para el registro de datos, utilizando el metodo de maxima verosimilitud (MLE).

Fit = markovchainFit(data = PE,confidencelevel = 0.95)
Fit

## $estimate
## MLE Fit 
##  A  3 - dimensional discrete Markov Chain defined by the following states: 
##  0, 1-5, 6+ 
##  The transition matrix  (by rows)  is defined as follows: 
##             0       1-5        6+
## 0   0.6605839 0.2299270 0.1094891
## 1-5 0.4625850 0.3061224 0.2312925
## 6+  0.1976285 0.3122530 0.4901186
## 
## 
## $standardError
##              0        1-5         6+
## 0   0.03471952 0.02048353 0.01413498
## 1-5 0.03966634 0.03226814 0.02804834
## 6+  0.02794888 0.03513120 0.04401395
## 
## $confidenceLevel
## [1] 0.95
## 
## $lowerEndpointMatrix
##             0       1-5        6+
## 0   0.5925349 0.1897800 0.0817850
## 1-5 0.3848404 0.2428780 0.1763188
## 6+  0.1428496 0.2433971 0.4038528
## 
## $upperEndpointMatrix
##             0       1-5        6+
## 0   0.7286330 0.2700740 0.1371931
## 1-5 0.5403296 0.3693669 0.2862663
## 6+  0.2524073 0.3811089 0.5763843
## 
## $logLikelihood
## [1] -1040.419

Dicha función devuelve una lista con todos los resultados de la estimaciones, incluyendo un objeto markovchain que posee la matriz de transición. Para el ejemplo anterior se tienen los siguientes resultados.

mc = Fit$estimate
summary(mc)

## MLE Fit  Markov chain that is composed by: 
## Closed classes: 
## 0 1-5 6+ 
## Recurrent classes: 
## {0,1-5,6+}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE

El diagrama de transición de la cadena es:

plot(mc)

Predicciones

Una de las utilidades más grandes de los Procesos Estocásticos, es su capacidad de hacer predicciones, esto se puede realizar mediante la función “predict()”, para eso observemos los ultimos registros del proceso

tail(PE)

## [1] "0"   "1-5" "0"   "6+"  "6+"  "1-5"

Podemos notar que la ultima semana llovió entre 1 y 5 veces, como el modelo es una cadena de markov, para realizar una predicción de la siguiente semana, necesita conocer la seana actual, las predicciones para las siguientes n=3 semanas dado que en la ultima vez se registraron entre 1 y 5 lluvias, son:

predict(mc,newdata="1-5",n.ahead=3)

## [1] "0" "0" "0"

U1A12 Estudio sobre las Cadenas de Markov en R

Isai Perez

2/3/2021