En este tutorial aplicaremos las técnicas para el análisis de tiempos promedios en cadenas de Markov. En primer lugar obtendremos los tiempos promedios acumulados de ocupación de los estados de una cadena de Markov en un horizonte de tiempo finito y después los tiempos de primera pasada. Finalmente, se presenta el análisis de tiempos promedio de absorción, lo cual es únicamente relevante para cadenas absorbentes.

Caso de aplicación: Modelo de mantenimiento de vías de trenes

En primer lugar, realizaremos el Problema 1 del archivo Complementaria 7.pdf que se encuentra en Sicua Plus.

Se modelará la situación como una cadena de Markov de tiempo discreto. Para esta situación se definen dos variables de estado:

\[\begin{eqnarray*} & Y_n = \text{ resultado de la revisión }n \\ & Z_n = \text{ resultado de la revisión }n-1\\ & X_n = (Z_n, Y_n) & \\ & S_X = \{(0,0),(0,1),(1,1),(1,2),(0,2),(2,0),(2,2)\} & \end{eqnarray*}\]

Donde \(0\) indica perfecto estado, \(1\) indica un daño de nivel amarillo y \(2\) un daño de nivel rojo.

Tenemos la siguiente matriz de probabilidades:

\[\textbf{P}= \begin{bmatrix} 0.6 & 0.25 & 0 & 0 & 0.15 & 0 & 0 \cr 0 & 0 & 0.7 & 0.3 & 0 & 0 & 0 \cr 0 & 0 & 0.4 & 0.6 & 0 & 0 & 0 \cr 0 & 0 & 0 & 0 & 0 & 0.5 & 0.5 \cr 0 & 0 & 0 & 0 & 0 & 0.5 & 0.5 \cr 1 & 0 & 0 & 0 & 0 & 0 & 0\cr 0 & 0 & 0 & 0 & 0 & 0.5 & 0.5\cr \end{bmatrix}\]

Vamos a crear el espacio de estados y la matriz de tasas de transición en R:

library(markovchain)
## Package:  markovchain
## Version:  0.6.9.12
## Date:     2018-08-23
## BugReport: http://github.com/spedygiorgio/markovchain/issues
estadosEjercicio = c("(0,0)","(0,1)","(1,1)","(1,2)","(0,2)","(2,0)","(2,2)")

probabilidadesEjercicio= matrix(data = c(0.6, 0.25,0,0,0.15,0,0,0,0,0.7,0.3,0,0,0,0,0,0.4,0.6,0,0,0,0,0,0,0,0,0.5,0.5,0,0,0,0,0,0.5,0.5,1,0,0,0,0,0,0,0,0,0,0,0,0.5,0.5), byrow = TRUE, nrow = 7, ncol=7)

probabilidadesEjercicio
##      [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,]  0.6 0.25  0.0  0.0 0.15  0.0  0.0
## [2,]  0.0 0.00  0.7  0.3 0.00  0.0  0.0
## [3,]  0.0 0.00  0.4  0.6 0.00  0.0  0.0
## [4,]  0.0 0.00  0.0  0.0 0.00  0.5  0.5
## [5,]  0.0 0.00  0.0  0.0 0.00  0.5  0.5
## [6,]  1.0 0.00  0.0  0.0 0.00  0.0  0.0
## [7,]  0.0 0.00  0.0  0.0 0.00  0.5  0.5

Ahora, creamos la cadena:

cmEjercicio = new(Class="markovchain", states = estadosEjercicio, byrow = TRUE,
transitionMatrix = probabilidadesEjercicio)

Vamos a verificar que la cadena sea ergódica. Para esto, hallamos el periodo y verificamos si es irreducible:

is.irreducible(cmEjercicio) 
## [1] TRUE
period(cmEjercicio)
## [1] 1

Ya que la cadena es ergódica, podemos hallar las probabilidades en estado estable:

steadyStates(cmEjercicio) 
##          (0,0)      (0,1)    (1,1)      (1,2)      (0,2)     (2,0)
## [1,] 0.3647416 0.09118541 0.106383 0.09118541 0.05471125 0.1458967
##          (2,2)
## [1,] 0.1458967

Análisis de tiempos promedio

En esta parte vamos a realizar analisis de tiempos de primera pasada y matriz de ocupación.

Cálculo de la matriz de ocupación

Si queremos conocer el número promedio de veces que la vía va a estar en perfecto estado en 15 revisiones si sabemos que hoy está en perfecto estado, ¿cómo lo podemos calcular?

La matriz de ocupación corresponde al número promedio que la cadena ocupa o visita el estado que se está analizando. Para obtener la matriz de ocupación en una CMTD elevamos M^t y sumamos la matriz identidad. Para el calculo vamos a definir una función que recibe como parámetro la matriz P y el tiempo n.

    matrizOcupacion <- function(mat, n){
      res = diag(rep(1, nrow(mat)))
      m = mat
      for (i in 1:n) {
        res = res + m
        m = m %*% m
      }
      return(res)
    }

matrizOcupacion(probabilidadesEjercicio, 15)

A partir de la matriz de ocupación sabemos que el número promedio de revisiones que las vía estará en perfecto estado dado que hoy estaba en perfecto estado es 6.6415.

Cálculo de tiempos de primera pasada

Si quisieramos saber cuántas revisiones pasarán hasta que la vía cumpla 2 revisiones en daños nivel rojo dado que ahora se encuentra en perfecto estado, ¿cómo lo podríamos calcular?

Debemos hallar el tiempo de primera pasada desde el estado (0,0) hasta el estado (2,2), es decir: \[ m_{(0,0),(2,2)}\] Podemos hallarlo utilizando la siguiente fórmula: \[m_{i,j} = p_{i,j} \sum\limits_{k\neq j} (1+m_{k,j})p_{i,k} \]

Otra opción para conocer el tiempo de primera pasada desde el estado (0,0) hasta el estado (2,2) es volver el estado (2,2) absorbente y calcular el tiempo antes de absorción para este estado. Esta operación es útil cuando el sistema de ecuaciones para hallar el timepo de primera pasada es muy grande.

En primer lugar obtenemos la matriz Q; submatriz de las probabilidades de transición entre estados transitorios, eliminado la fila y columna relacionada con el estado absorbente que vamos a crear. Luego hallamos la matriz B; número promedio de visitas antes de absorción, restando a la matriz identidad I la matriz Q y hallando la inversa de la matriz resultante. Por último, sumamos el número promedio de semanas que la cadena va a estar en cada uno de los estados transitorios iniciando en el estado (0,0) antes de ser absorbida por el estado (2,2).

#Posición del estado que va a ser absorbente
indiceAbsorbente=7
#Matriz Q
matrizQ<-probabilidadesEjercicio[-indiceAbsorbente,-indiceAbsorbente]
#Matriz identidad
matrizI<-diag(nrow(matrizQ))
#Matriz B
matrizB<-solve(matrizI-matrizQ)

revisionesAntes22<-sum(matrizB[1,])

revisionesAntes22

A pertir de los cálculos anteriores podemos observar que el número de semanas promedio que pasarán desde que la vía complete 2 semanas en daños de nivel rojo dado que lleva 2 semanas en perfecto estado es de 10.7083 semanas.