El objetivo de este tutorial es que aprendan a definir y manejar matrices en R. El uso de matrices aplicado a la construcción de Cadenas de Markov también será tratado a lo largo del tutorial.

Definición de matrices

Una matriz es un arreglo de datos de dimensiones mxn, en donde m representa la cantidad de filas y n la cantidad de columnas. En general, una matriz en R es la unión de una serie de vectores del mismo tamaño. Se utiliza de forma similar a un Data Frame visto en la semana 3. Utilizamos la función matrix() en donde incluímos los valores de los datos que van a conformar nuestra matriz y el número de filas y columnas que tendrá. Debemos tener en cuenta que la matriz se crea por columnas por default a menos que lo especifiquemos con el argumento byrow.

¡Investiga!
¿Cuál es la diferencia entre una matriz y un data frame en R?

mat = matrix(c(1,4,1,3,2,4,4,0,3,2,1,1,4,2,2,0),nrow=4,ncol=4,byrow = TRUE)

También podemos asignarle nombres a las columnas y filas de las matrices. Además se puede utilizar el comando colnames() para asignar nombre sólo a las columnas y rownames() para asignar nombre a las filas.

dimnames(mat) = list(c("SixPack Team","Alianza Team","Real FC","Fortaleza FC"),c("Pos","Goles Marcados","Goles Recibidos","Diferencia"))
##              Pos Goles Marcados Goles Recibidos Diferencia
## SixPack Team   1              4               1          3
## Alianza Team   2              4               4          0
## Real FC        3              2               1          1
## Fortaleza FC   4              2               2          0

Otra forma de crear matrices es uniendo una serie de vectores del mismo tamaño. Estos pueden unirse por columnas cbind() o por filas rbind():

a = seq_len(4)
b = rep(c(4,2),each = 2)
c = c(1,4,1,2)
d = b - c
mat1 = cbind(a,b,c,d)
##      a b c d
## [1,] 1 4 1 3
## [2,] 2 4 4 0
## [3,] 3 2 1 1
## [4,] 4 2 2 0

Otro ejemplo, es si queremos adicionar un nuevo valor a nuestra matriz inicial mat:

Mullaneiros = c(5,0,3,-3)
mat = rbind(mat,Mullaneiros)
##              Pos Goles Marcados Goles Recibidos Diferencia
## SixPack Team   1              4               1          3
## Alianza Team   2              4               4          0
## Real FC        3              2               1          1
## Fortaleza FC   4              2               2          0
## Mullaneiros    5              0               3         -3

En general, se accede a los elementos de una matriz de la misma forma en que se hace con un Data Frame. En este caso se puede poner el nombre de la fila/columna o el índice que hace referencia a esta:

mat["Fortaleza FC",]
##             Pos  Goles Marcados Goles Recibidos      Diferencia 
##               4               2               2               0
mat[,4]
## SixPack Team Alianza Team      Real FC Fortaleza FC  Mullaneiros 
##            3            0            1            0           -3

¡Ponte a prueba!
¿Cómo podrías obtener la posición y los goles marcados únicamente para Fortaleza FC y Real FC?

Una función muy útil en el manejo de matrices es apply(), que aplica otra función sobre los valores de una matriz. En esta función el primer argumento es el conjunto de datos sobre el que se quiere aplicar la función, el segundo es 1 o 2 si se hace sobre filas o columnas, respectivamente; y el tercero, es la función que se quiere aplicar. Por ejemplo, si queremos calcular la media de las columnas de la matriz mat:

apply(mat, 2, mean)
##             Pos  Goles Marcados Goles Recibidos      Diferencia 
##             3.0             2.4             2.2             0.2

¡Investiga!
¿Cuáles otras funciones puedo aplicar en una matriz con la función “apply()”?

Para evaluar si un conjunto de datos es de tipo matriz se utiliza el comando is.matrix() . Con as.matrix() podemos convertir un arreglo de datos a un objeto que sea de tipo matriz. Esto es importante debido a que en R se pueden hacer operaciones particulares con las matrices.

Operaciones entre matrices

En este caso vamos a crear dos matrices llenas de números aleatorios generados por el comando rnorm().

A = matrix(rnorm(n = 36,mean = 20,sd = 8),nrow=6,ncol=6)
B = matrix(rnorm(n = 36,mean = 10,sd = 3),nrow=6,ncol=6)

Las operaciones más comunes de álgebra lineal y manejo de matrices se pueden evaluar con los siguientes comandos:

Comando Operación
t(A) Transpuesta de la matriz A
A + B Suma de matrices A y B
A - B Resta de las matrices A y B
A%*%B Producto matricial entre A y B
A*B Matriz producto elemento por elemento de A y B
det(A) Determinante de la matriz A
solve(A) Inversa de la matriz A
solve(A,b) Solución del sistema de ecuaciones Ax = b
eigen(A) Valores y vectores propios de A
diag(A) Diagonal de la matriz A
colSums(A) Suma de cada columna de A
rowSums(A) Suma de cada fila de A

Vamos a utilizar algunas de las funciones más comunes sobre las matrices creadas:

colSums(A%*%B)
## [1] 9313.894 6611.168 7736.774 8186.777 9745.369 6725.948
rowSums(A+B)
## [1] 162.2711 148.7829 186.6081 226.8588 183.6227 212.8624

Cadenas de Markov

El uso de matrices nos sirve para definir la representación de una Cadena de Markov en R. Por ejemplo, en las cadenas continuas podemos definir la matriz generadora (o de tasas de transición) Q y en las discretas definimos la matriz de probabilidades de transición a un paso P. En R se pueden encontrar muchos paquetes que permiten manejar y desarrollar Cadenas de Markov, entre estos están “cmtcmove”, “DOBAD”, “markovchain” y “DTMCPack”.

Para aprender a solucionar problemas utilizando cadenas de Markov utilizaremos el Problema 1 del archivo “Complementaria 3 (Q).pdf” que está en SICUA+.

¡Ponte a prueba!
Modela el sistema de inventario de Cabbage Patch Kids como una cadena de Markov.

Luego de hacer la modelación, pasamos a R para resolver el problema con ayuda de herramientas computacionales. Primero debemos definir la matriz P (CMTD). Para estimar las probabilidades de un Proceso de Poisson utilizamos la función dpois() que tiene como parámetro el valor de x y el valor de \(\lambda\).

\[P(N(t) = x)=\frac{e^{-\lambda t} \cdot (\lambda t)^{x}}{x!}\]

Si queremos calcular el valor de la probabilidad acumulada:

\[P(N(t) \leq x)=\sum_{i=0}^{x}\frac{e^{-\lambda t} \cdot (\lambda t)^{i}}{i!}\]

utilizamos la función ppois() que tiene como parámetro el valor de x y el valor de \(\lambda\), también, se puede configurar para que se devuelva el valor de la cola derecha \(P(N(t) > x)\) utilizando el parámetro lower.tail. Esto se podrá de la misma forma con una suma de las probabilidades individuales calculadas con dpois(). Probamos con el valor calculado por ambas funciones y comparamos con los valores calculados manualmente:

a = exp(-3)*(3^0)/factorial(0)
b = exp(-3)*(3^1)/factorial(1)

Calculamos los valores de \(P(N(1)=0)\), \(P(N(1)=1)\) y \((P(N(1) \leq 1)=(P(N(1)=0)+P(N(1)=1)))\), que en este caso corresponden a a, b y a+b, respectivamente:

## [1] 0.04978707
## [1] 0.1493612
## [1] 0.1991483

Ahora, comparamos con los valores calculados por las dos funciones a utilizar:

dpois(0,lambda=3)
## [1] 0.04978707
dpois(1,lambda=3)
## [1] 0.1493612
ppois(1,lambda=3)
## [1] 0.1991483
ppois(1,lambda=3,lower.tail = FALSE)
## [1] 0.8008517

¡Ponte a prueba!
Ya que sabemos cómo utilizar las funciones que estiman las probabilidades de un Proceso de Poisson, crea la matriz P de nuestra CMTD.

Para conocer el número de muñecos que habrá en inventario el viernes dentro de 4 semanas, elevamos la matriz P a la cuarta potencia y calculamos el valor deseado:

p2 = p %*% p
p4 = p2 %*% p2
p4[3,1]*2 + p4[3,2]*3 + p4[3,3]*4 + p4[3,4]*5
## [1] 2.761991

¡Ponte a prueba!
¿Cómo podrías hallar las probabilidades en estado estable para el nivel de inventario de Cabbage Patch Kids?

Paquete markovchain

El paquete markovchain proporciona un conjunto de métodos que permiten manejar cadenas de Markov de tiempo discreto y contínuo.

Cadenas de Markov de tiempo discreto

Para utilizar este paquete se debe llamar e instalar, luego es necesario crear un objeto de la clase markovchain, el cual está compuesto por los estados y la matriz de transición de probabilidad.

library(markovchain)
## Package:  markovchain
## Version:  0.6.9.8-1
## Date:     2017-08-15
## BugReport: http://github.com/spedygiorgio/markovchain/issues
estados<-c("inv0","inv1","inv2","inv3","inv4","inv5")
cmtdCabbage<-new("markovchain", states=estados, transitionMatrix=p)
print(cmtdCabbage)
##           inv0      inv1       inv2       inv3       inv4       inv5
## inv0 0.1847368 0.1680314 0.22404181 0.22404181 0.14936121 0.04978707
## inv1 0.1847368 0.1680314 0.22404181 0.22404181 0.14936121 0.04978707
## inv2 0.8008517 0.1493612 0.04978707 0.00000000 0.00000000 0.00000000
## inv3 0.5768099 0.2240418 0.14936121 0.04978707 0.00000000 0.00000000
## inv4 0.3527681 0.2240418 0.22404181 0.14936121 0.04978707 0.00000000
## inv5 0.1847368 0.1680314 0.22404181 0.22404181 0.14936121 0.04978707

Para visualizar la cadena de Markov se puede utilizar la función “plot()”.

plot(cmtdCabbage)

Por último, el paquete markovchain ofrece múltiples funciones que se pueden utilizar para analizar cadenas de Markov, por ejemplo; se puede saber si una cadena es irreducible, se puede obtener el vector de probabilidades en estado estable, entre otras.

¡Investiga!
¿Cuáles funciones ofrece el paquete markovchain para analizar cadenas de Markov?

Cadenas de Markov de tiempo contínuo

La creación de cadenas de Markov de tiempo contínuo es similar al proceso anterior, después de instalar e importar el paquete markovchain se debe crear la cmtc. Para ilustrar utilizaremos como ejemplo el funcionamiento de un ascensor. Este ascensor puede estar funcionando o dañado en cada momento del tiempo. A continuación creamos la matriz de tasas de transición para modelar los estados del ascensor.

estadosAscensor <- c("Funcionando", "Dañado")
byRow <- TRUE
matrizQAscensor <- matrix(data = c(-3, 3,1, -1), nrow = 2, byrow = byRow, dimnames = list(estadosAscensor, estadosAscensor))

Para crear el objeto de tipo cadena de Markov utilizamos la función “new()” pero ahora el primer parámetro es “ctmc”, este indica que estamos creando una cadena de Markov de tiempo contínuo.

CTMCascensor <- new("ctmc", states = estadosAscensor,byrow = byRow, generator = matrizQAscensor,
name = "Matriz de transición ascensor")

¡Ponte a prueba!
Formula e implementa el Problema 2 del archivo“Complementaria 3 (Q).pdf” que está en SICUA+