Instalacion del paquete.

Primero se debe cargar e instalar el paquete markovchain desde CRAN. Para hacerlo se usa la función install.packages(arg) como se muestra a continuacion:

install.packages('markovchain')

Alternativamente, si se trabaja con una IDE, como RStudio o Pycharm, es posible instalar el módulo usando el menu de paquetes de la aplicación. En el caso de RStudio, el IDE más popular para R, se deben completar los siguientes pasos

1. Dar click en Tools → Install Packages
2. Seleccionar Repository (CRAN) en el espacio Install from:
3. Ingresar el nombre del paquete (o varios nombres, separados por un espacio o coma), en este caso markovchain
4. Dejar la opción Install dependencies seleccionada por default.
5. Dar click en Install

Finalmente, cargamos el módulo usando la función library() o require(), como se muestra a continuación. Esto se puede hacer en la consola o directamente en el script o rmd.

library(markovchain)

Despues de correr el codigo, lo siquiente debería aparecer en la consola. En caso contrario, el paquete no fue instalado correctamente.

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

En general se recomienda más el uso de library() sobre require(). Para más información sobre la diferencia entre estas dos funciones revisar este link.

Uso del paquete.

A continuación se muestra como hacer uso del paquete para crear un objeto markovchain, obtener un resumen de sus principales características, clasificar los estados, encontrar probabilidades de transición, obtener la matriz de trancisión en \(n\) pasos, obtener la distribución de la variable \(X_n\), calcular una distribución estacionaria y obtener el gráfo de la cadena. Para más información sobre el paquete revisar este link.

Para ejemplificar, se supondrá una cadena de markov \(X = \{X_1, X_2, ...\}\) definida en un espacio de estados \(S = \{a, b, c\}\) y cuya matriz de trancisión esta dada por: \[P = \begin{bmatrix} 0 & 0.5 & 0.5\\ 0.5 & 0 & 0.5\\ 0.5 & 0.5 & 0 \end{bmatrix}\]

Creación de una CM.

La primera forma de crear una instancia de markovchain es definir una matriz de trancisión P y un vector de estados S, esto se puede hacer usando las funciones matrix() y letters respectivamente como se muestra a continuación.

P <- matrix(c(0, 0.5, 0.5,
              0.5, 0, 0.5,
              0.5, 0.5, 0), nrow = 3, byrow = T);
S <- letters[1:3];

Aquí la función matrix() toma tres argumentos: El primero es un vector que contiene las entradas de la matriz, el segundo nrow es el numero deseado de columnas y byrow indica si la matriz se va a llenar por columnas (FALSE) o por filas (TRUE). En el caso de S la función letters nos da un vector con las letras empezando por la primera (a) hasta la n-esima letra.

Podemos usar la función show() para mostrar como se ve nuestra matriz y el conjunto de estados.

show(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

show(S)

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

Ahora podemos crear una instancia de markovchain usando la función new(). La clase markovchain tiene 4 atributos: transitionMatrix (la matriz de trancisión), states (los estados), name (el nombre del objeto) y byrows (Si la matriz es estocástica por filas (TRUE) o columnas (FALSE)). A continuación creamos la cadena de markov del ejemplo.

CM_1 <- new("markovchain", transitionMatrix = P, states = S, name = "Ejemplo 1", byrow = TRUE);
show(CM_1)

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

En este caso show() muestra la “forma” del objeto markovchain o bien, su definición como tipo str.

P_alt1 <- matrix(c(0, 0.5, 0.5,
              0.5, 0, 0.5,
              0.5, 0.5, 0), nrow = 3, byrow = T);
colnames(P_alt1) <- rownames(P_alt1) <- letters[1:3];
CM_alt1 <- as(P_alt1, "markovchain")
show(CM_alt1)

## Unnamed Markov chain 
##  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

P_alt2 <- matrix(c(0, 0.5, 0.5,
              0.5, 0, 0.5,
              0.5, 0.5, 0),
                 nrow = 3,
                 byrow = T,
                 dimnames = list(letters[1:3], letters[1:3]))

CM_alt2 <- new("markovchain", transitionMatrix = P_alt2, name = "Instancia Alternativa", byrow = TRUE);
show(CM_alt2)

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

Características principales de la CM.

La función summary() da un resumen de las principales características de la cadena de Markov. En este resumen podemos ver las clases cerradas, recurrentes y transitivas, asi como si la CM es irreducible y cuales son sus estados absorbentes.

summary(CM_1)

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

De igual forma, podemos solicitar las distintas partes de este resumen por separado usando las funciones recurrentClasses(), transientStates(), absorbingStates(), is.irreducible(), communicatingClasses().

• Clases recurrentes

recurrentClasses(CM_1)

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

• Estados transitivos

transientClasses(CM_1)

## list()

• Estados absorbentes

absorbingStates(CM_1)

## character(0)

• Clases comunicativas

communicatingClasses(CM_1)

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

• Si la CM es irreducible

is.irreducible(CM_1)

## [1] TRUE

• El periodo

period(CM_1)

## [1] 1

Tambien podemos revisar si un estado es accesible desde otro usando is.accesible(object, from, to) como se muestra a continuación.

is.accessible(CM_1, "a", "c")

## [1] TRUE

Gráfo de la CM

Para visualizar el diagrama de transición de la CM se puede usar la función plot(). El gráfo de la CM es un objeto tipo igraph asi que puede recibir como argumentos adicionales a la CM cualquier argumento que recibiría un gráfico de este tipo. Para más información sobre los parámetros que se pueden especificar revisar el manual del paquete igraph.

plot(CM_1,
    edge.arrow.size = 0.75,
    vertex.size = 20, edge.label.cex = 2,
    vertex.label.cex = 2)

Análisis probabilístico

Para conocer la probabilidad de trancisión en 1 paso entre un estado y otro basta con utilizar la función transitionProbability(object, t0, t1), en donde object es la CM, t0 es el estado en el tiempo 0 y t1 el estado en el tiempo 1.

Por ejemplo, la probabilidad de transición en un paso del estado \(a\) al estado \(b\) sería:

transitionProbability(CM_1, 'a', 'c')

## [1] 0.5

El resultado también se puede obtener dado que se sabe que eso corresponde a la entrada \(P[1, 3]\) de la patriz \(P\).

show(CM_1[1, 3])

## [1] 0.5

Así mismo podemos convertir la markovchain a un dataframe que contenga las probabilidades de trancisión en un paso de cada estado a otro estado.

mcDf <- as(CM_1, "data.frame");
show(mcDf)

##   t0 t1 prob
## 1  a  a  0.0
## 2  a  b  0.5
## 3  a  c  0.5
## 4  b  a  0.5
## 5  b  b  0.0
## 6  b  c  0.5
## 7  c  a  0.5
## 8  c  b  0.5
## 9  c  c  0.0

CM_5 <- CM_1^5;
show(CM_5)

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

transitionProbability(CM_5, "a", "c")

## [1] 0.34375

show(CM_5[1, 3])

## [1] 0.34375

mcDf5 <- as(CM_5, "data.frame");
show(mcDf5)

##   t0 t1    prob
## 1  a  a 0.31250
## 2  a  b 0.34375
## 3  a  c 0.34375
## 4  b  a 0.34375
## 5  b  b 0.31250
## 6  b  c 0.34375
## 7  c  a 0.34375
## 8  c  b 0.34375
## 9  c  c 0.31250

Distribución de la variable \(X_n\)

Usando las ecuaciones de Chapman-Kolmogorov se sabe que \[\pi_n = \pi_0P_n = \pi_0P^n\] Así, si el vector inicial de la cadena de Markov CM_1 esta dado por \[\pi_0 = (0.5, 0.2, 0.3)\] entonces podemos calcular \(\pi_6\), por ejemplo, como se muestra a continuación

pi_0 <- c(0.5, 0.2, 0.3);
pi_6 <- pi_0*CM_1^6;
show(pi_6)

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

Y podemos confirmar que en efecto el vector tiene entradas positivas y suma 1.

show(sum(pi_6))

## [1] 1

show(all(pi_6 >= 0))

## [1] TRUE

Para facilitar este tipo de análisis en un futuro, es conveniente hacer una función que tome como parametros una distribucion inicial \(\pi_0\) y una \(n\). A continuación se muestra una implementación de dicha función.

pi <- function(CM, pi_0, n = 1) {
    return(pi_0 * CM^n)
}

Asi pues, \(\pi_6\) sería

pi(CM_1, pi_0, 6)

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

Distribución estacionaria

Una distribución es estacionaria si sucede que \[\pi^{*} = \pi^{*}P\] Se puede encontrar una distribución estacionaria para una CM utilizando la función steadyStates() como se muestra a continuación.

steadyStates(CM_1)

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

Ejercicio 1

Para el primer ejercicio se pide determinar las clases comunicantes, diagrama de comunicación, los estados cerrados, recurrentes y transitivos y el periodo de varias cadenas de markov dadas por una matriz de trancisión distinta.

En este caso las funciones que se utilizarán son: - communicatingClasses() para obtener las clases comunicantes. - summary() para obtener los estados cerrados, recurrentes y transitivos. - period() para obtener el periodo. - plot() para obtener el diagrama.

Para no estar repitiendo los mismos comandos una y otra vez, conviene definir una función que nos permita obtener toda esta información de las distintas CM.

ejercicio1 <- function(P){
    cm <- as(P, "markovchain")
    print("The communicating classes are:")
    show(communicatingClasses(cm))
    summary(cm)
    if (is.irreducible(cm)){
        print("The period is:")
        show(period(cm))
    }

    plot(cm, edge.arrow.size = 0.75,
        vertex.size = 20, edge.label.cex = 1,
        vertex.label.cex = 1)
}

A continuación se muestran los ejercicios y su solución.

P1 <- matrix(c(1, 0, 0, 0,
                0.5, 0, 0.5, 0,
                0, 0.3, 0.7, 0,
                0.2, 0, 0.8, 0),
                nrow = 4,
                byrow = T,
                dimnames = list(letters[1:4], letters[1:4]))
ejercicio1(P1)

## [1] "The communicating classes are:"
## [[1]]
## [1] "a"
## 
## [[2]]
## [1] "b" "c"
## 
## [[3]]
## [1] "d"
## 
## Unnamed Markov chain  Markov chain that is composed by: 
## Closed classes: 
## a 
## Recurrent classes: 
## {a}
## Transient classes: 
## {b,c},{d}
## The Markov chain is not irreducible 
## The absorbing states are: a

P2 <- matrix(c(0.3, 0, 0.7,
                0.4, 0.6, 0,
                0.5, 0, 0.5),
                nrow = 3,
                byrow = T,
                dimnames = list(letters[1:3], letters[1:3]))
ejercicio1(P2)

## [1] "The communicating classes are:"
## [[1]]
## [1] "a" "c"
## 
## [[2]]
## [1] "b"
## 
## Unnamed Markov chain  Markov chain that is composed by: 
## Closed classes: 
## a c 
## Recurrent classes: 
## {a,c}
## Transient classes: 
## {b}
## The Markov chain is not irreducible 
## The absorbing states are: NONE

P3 <- matrix(c(0, 0.2, 0, 0.8,
                0.4, 0.6, 0, 0,
                0, 0, 1, 0,
                0, 0, 0.3, 0.7),
                nrow = 4,
                byrow = T,
                dimnames = list(letters[1:4], letters[1:4]))
ejercicio1(P3)

## [1] "The communicating classes are:"
## [[1]]
## [1] "a" "b"
## 
## [[2]]
## [1] "c"
## 
## [[3]]
## [1] "d"
## 
## Unnamed Markov chain  Markov chain that is composed by: 
## Closed classes: 
## c 
## Recurrent classes: 
## {c}
## Transient classes: 
## {a,b},{d}
## The Markov chain is not irreducible 
## The absorbing states are: c

P4 <- matrix(c(0.25, 0.25, 0.25, 0.25,
                0, 0.5, 0.5, 0,
                0, 0, 0, 1,
                0, 0, 1, 0),
                nrow = 4,
                byrow = T,
                dimnames = list(letters[1:4], letters[1:4]))

ejercicio1(P4)

## [1] "The communicating classes are:"
## [[1]]
## [1] "a"
## 
## [[2]]
## [1] "b"
## 
## [[3]]
## [1] "c" "d"
## 
## Unnamed Markov chain  Markov chain that is composed by: 
## Closed classes: 
## c d 
## Recurrent classes: 
## {c,d}
## Transient classes: 
## {a},{b}
## The Markov chain is not irreducible 
## The absorbing states are: NONE

P5 <- matrix(c(0, 0, 0, 0, 1,
                0, 0, 1, 0, 0,
                0, 0, 0, 0.5, 0.5,
                0, 0, 1, 0, 0,
                0.2, 0.2, 0.2, 0.2, 0.2),
                nrow = 5,
                byrow = T,
                dimnames = list(letters[1:5], letters[1:5]))

ejercicio1(P5)

## [1] "The communicating classes are:"
## [[1]]
## [1] "a" "b" "c" "d" "e"
## 
## Unnamed Markov chain  Markov chain that is composed by: 
## Closed classes: 
## a b c d e 
## Recurrent classes: 
## {a,b,c,d,e}
## Transient classes: 
## NONE 
## The Markov chain is irreducible 
## The absorbing states are: NONE
## [1] "The period is:"
## [1] 1

Ejercicio 2

Suponga que un mexicano puede estar en 1 de 4 estados: Rico, Clase media, Pobre o Pobreza extrema. Suponga las siquientes probabilidades de transición:

• Si un mexicano es Rico, su siguiente generación será:
  • Rico: \(0.0\)
  • Clase media: \(0.75\)
  • Pobre: \(0.2\)
  • Pobreza extrema: \(0.05\)
• Si un mexicano es de Clase media, su siguiente generación será:
  • Rico: \(0.05\)
  • Clase media: \(0.2\)
  • Pobre: \(0.3\)
  • Pobreza extrema: \(0.45\)
• Si un mexicano es Pobre, su siguiente generación será:
  • Rico: \(0.1\)
  • Clase media: \(0.4\)
  • Pobre: \(0.3\)
  • Pobreza extrema: \(0.2\)
• Si un mexicano está en Pobreza extrema, su siguiente generación será:
  • Rico: \(0.27\)
  • Clase media: \(0.15\)
  • Pobre: \(0.03\)
  • Pobreza extrema: \(0.55\)

S <- c("Rico", "Clase Media", "Pobre", "Pobreza extrema");
P <- matrix(c(0, 0.75, 0.2, 0.05,
            0.05, 0.2, 0.3, 0.45,
            0.1, 0.4, 0.3, 0.2,
            0.27, 0.15, 0.03, 0.55),
            nrow = 4,
            byrow = T,
            dimnames = list(S, S));

CM <- new("markovchain", transitionMatrix = P, name = "CM de Riqueza en México", byrow = T);

transitionProbability(CM, "Clase Media", "Rico")

## [1] 0.05

transitionProbability(CM^2, "Clase Media", "Rico")

## [1] 0.1615

transitionProbability(CM^3, "Clase Media", "Rico")

## [1] 0.1386

transitionProbability(CM, "Pobreza extrema", "Clase Media")

## [1] 0.15

transitionProbability(CM^2, "Pobreza extrema", "Clase Media")

## [1] 0.327

transitionProbability(CM^3, "Pobreza extrema", "Clase Media")

## [1] 0.292875

steadyStates(CM)

##           Rico Clase Media     Pobre Pobreza extrema
## [1,] 0.1376442   0.2925528 0.1813545       0.3884486