Cadenas de Markov (Markov chains)

Cadenas de Markov implementadas en R

• Instalar el paquete markovchain y activarlo

library(markovchain)

## Package:  markovchain
## Version:  0.8.5-2
## Date:     2020-09-07
## BugReport: https://github.com/spedygiorgio/markovchain/issues

• La libreria nos proveera objetos para realizar análisis estadísticos de cadenas de markov a tiempos discretos. Asumamos que tenemos una cadena de markov x={x1,x2...} definida en el espacio de estado S={a,b,c} y cuya transición es:

\[ P = \left( {\being{array}{ccc} 0 & .5 & .5 \\ .5 & 0 & .5 \\ .5 & .5 & 0 \\ \end{array} } \right) \]

1.- Crear matriz de transición P:

P= matrix(c(0,.5,.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

#nrow: es para que ponga el número de filas
#byrow: es para que se almacenen los datos de la matriz
#c(): es para concatenar

• Crear la matriz de transición, creamos objeto "markovchain" así:

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

## Cadena 1 
##  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.0 0.5 0.5
## b 0.5 0.0 0.5
## c 0.5 0.5 0.0

• La estructura del objeto "mc" es:

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"

• Sumario de 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

• Visualizar de manera gráfica 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 obtenida

recurrentClasses(mc)

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

transientStates(mc)

## character(0)

absorbingStates(mc)

## character(0)

Ejercicio para predecir el clima

• Se presenta el caso de una ciudad desconocida donde llueve n número de veces a la semana, los datos estan 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 tabla de frecuencias para rain es:

table(rain$rain)

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

• Observar los primeros registros:

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

• Se propone de modelo una cadena de markov finita de 3 estados para visualizar las veces que llueve en la semana

PE <- rain$rain
head(PE)

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

• La funcion "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

• La función "markovchainFit()" estima la matriz de transición para el registro de datos, haciendo uso del método de máxima verosimilitud (MLE)

Fit <- markovchainFit(data=PE, confidencelevel = .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 de todos los resultados estimados. Para el ejemplo anterior se tienen los diguientes resultados:

mc <- Fit$estimate
mc

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

#Porcentaje de, al pasar un evento, suceda el otro tomando en cuenta el evento anterior inmediato

• Gráfica de la cadena de Markov

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 Se pueden realizar predicciones a tráves de la función "predict()". Con los últimos registros del proceso, ejemplificaremos esto

tail(PE)

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

• Se mira que la última semana llovio de 1 a 5 veces. Para realizar la predicción de la siguiente semana con la cadena de markov, es necesario tener el valor de la semana actual. Para las siguentes n=3 semanas, las predicciones tomando en cuenta la semana actual serían:

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

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