Cadenas de Markov

Cadenas de Markov en R

Proceso Estocásticos

Instalación del Paquete

Para instalar un paquete en R se utiliza el comando “install.packages()”, para mayor información se digita en el compilador el código:

help(“install.packages”)

De esa forma se accede a la documentacion de dicha funcion almacenada en el CRAN. Para analizar cadenas de markov finitas utilizaremos el paquete “markovchain” version 0.6.9.11 publicado por: Giorgio Alfredo Spedicato, Tae Seung Kang y otros. En el siguiente chunk se presenta los comandos necesarios para instalar y cargar al sistema el paquete. Note que la primera linea de código tiene el simbolo “#” que es utilizado para comentar.

# install.packages("markovchain")

library(markovchain)

## Package:  markovchain
## Version:  0.8.6
## Date:     2021-05-17
## BugReport: https://github.com/spedygiorgio/markovchain/issues

El comando “library” carga al sistema el paquete

Introducción a markovchain

Para conocer los detalles, especificaciones o ejemplos de la libreria basta con digitar el siguiente comando, que automaticamente los redireccionará a su propia documentación.

  help("markovchain")

Esta libreria pretende proveer objetos para realizar analisis estadísticos de cadenas de markov a tiempos discretos. Asumamos que tenemos una cadena de markov X={X1,X2,…} definida en el espacio de estados S={a,b,c} y cuya matriz de transición es:

\[P=\lbrace 00.50.50.500.50.50.50 \rbrace\]

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

Crear la matriz de transicion P:

P = matrix(c(0,0.5,0.5,.5,0,.5,.5,.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 “nrows” de la funcion matrix es para declarar el numero de filas que deseamos que nuestra matriz P posea, y el argumento “byrows” es para que almacene los elementos de la matriz almacenados en c(), fila por fila.

Crear la matriz de transición creamos el objeto “markovchain” de la siguiente forma: mc = new(“markovchain”,transitionMatrix=P,states=c(“a”,“b”,“c”),name=“Cadena 1”) Una revisión previa al análisis de nuestra cadena se puede realizar mediante los comandos “str()” y “summary”, que devuelven la estructura del objeto y el resumen general de los resultados respectivamente. Para mayor informacion revisar los comandos mediante la función help().

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

La estructura del objeto mc es:

library(markovchain)
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 tippo markovchain, que representa una lista con los argumentos provistos declarados en el constructor.

El resumen de la cadena es:

summary(mc)

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 visualizar el diagrama de transición de la cadena lo realizamos con la función plot(mc)

plot(mc)

Otras funciones importantes 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:

recurrentClasses(mc)

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

recurrentClasses(mc)

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

transientStates(mc)

## character(0)

absorbingStates(mc)

## character(0)

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 probabilidad de transicion en un paso del estado “a” al estado “c” es:

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

## [1] 0.5

Recuerde que dicha probabilidad es un elemento de la matriz de transición P, por lo tanto, la probabilidad de transicion del estado “a” al estado “b” es simplemente P23

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 # El número de pasos al futuro mc ^ n

n = 5 # El número de pasos al futuro
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 pueden conocer la distribución de la cadena en n pasos adelante (P(Xn)) multiplicando la distribucion inicial de X0 por la matriz de transición en n pasos (Pn), calcule la distribución de la cadena en el tiempo n = 6, si la ditribución inicial de la cadena es “(0.5, 0.2, 0.3)”.

X0 = c(0.5,0.2,0.3) # La distribucion de X en t = 0
n = 6
Xn = X0*(mc^n)
Xn

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

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

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 distribución 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

Recuerde que los tiempos medio de recurrencia son los inversos multiplicativos de la distribución estacionaria y pueden 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=⎛⎝⎜⎜⎜⎜1/2001/41/21/21/21/401/21/21/40001/4⎞⎠⎟⎟⎟⎟\] Para resolver este problema del segundo examen basta con almacenar correctamente la matriz de transición y usar los comandos presentados.

Q = matrix(data = 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)
Q

##      [,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

Por lo tanto la cadena de markov clasificada es:

Ejer1 = new("markovchain",transitionMatrix = Q,states = c("0","1","2","3"))
summary(Ejer1)

## 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(Ejer1)

Y la distribución estacionaria de la matríz es:

steadyStates(Ejer1)

##      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, (Inferencia y Estadistica 2).

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 de la forma, 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 transicion 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_{ij}\) 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)

Y la distribución estacionaria de la matriz es:

steadyStates(mc)

##              0       1-5        6+
## [1,] 0.5008871 0.2693656 0.2297473

##Predicciones

Una de las utilidades mas grandes de los Procesos Estocásticos, es su capacidad de realizar 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"

Se observa que la última 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 semana actual, las prediciones para las siguientes* n = 3* semanas dado que en la utlima vez se registraron entre 1 y 5 lluvias, son:

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

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

Para mayor información de la funcionalidad del paquete o mayor conocimiento de otras funciones pueden revisar el artículo de Giorgio Alfredo, Tae Seung y otros llamado “The markovchain Package: A Package for Easily Handling Discrete Markov Chains in R” mediante el siguiente Enlace