4.6 Probabilidades de transición

El aspecto funcdamental para estudiar el comportamiento de cualquier CMTC es la obtención y análisis de las probabilidades de transición entre los estados del procso a partir de la matriz de tasas obtenida en puntos anteriores. Aunque haremos un desarrollo teórico de este rporblema, veremos que es necesario un algoritmo de computación para obtener dichas probabilidades.

Sea \(X(t), t \geq 0\) una CMTC con espacio de estados \(S = \{1,2,...,N\}\) y con matriz de tasas \(R = [r_{ij}]\). Si asumimos que la distribución de probabilidad en el estado incial, \(X(0)\) es conocida, entonces tenemos que:

\[P(X(t) = j) = \sum_{i=1}^N P(X(t) = j | X_0 = i)P(X_0 = i), \quad 1 \leq j \leq N.\] Es necesario obtener \(P(X(t) = j | X_0 = i) = p_{ij}(t)\) para obtener la función de distribución de probabilidad de \(X(t)\). Antes de ver como obtener dichas probabildiades introducimos la notación necesaria.

ya hemos visto antes que una CMTC permance un tiempo \(Exp(r_i)\) en el estado \(i\) con \(r_i = \sum_{j=1}^N r_{ij}\) y, si \(r_i > 0\) entonces pasamos al estado \(j\) con probabilidad \(p_{ij} = r_{ij}/r_i\). Asumiendo que existe un número finito \(r\) que satisface:

\[r \geq max(r_i), \quad 1\leq i \leq N\] podemos definir la matriz \(\hat{P} = [\hat{p}_{ij}]\) como:

\[\begin{equation} \hat{p}_{ij} = \begin{cases} 1-\frac{r_i/r} & \text{ si } i=j\\ r_{ij}/r & \text{ si } i \neq j \end{cases} \end{equation}\]

que es una matriz estocástica.

Teorema 4.2 La matriz de tarnsición de probabilidades, \(P\), de una CMTC viene dada por:

\[P(t) = \sum_{k=0}^{\infty} e^{-rt}\frac{(rt)^k}{k!} \hat{P}^k\]

Esta forma de obtener \(P\) se denomina método de uniformización, y proporciona un método numérico para obtener las probabilidades de transición deseadas, sin más que usar los \(m\) primeros términos de la serie infinita definida. Para asegurar la convergencia de dichas probabilidades se usan como reglas habituales:

\[m \approx max\{rt + r\sqrt{rt}, 20\},\]

y

\[r = max(r_i), \quad 1\leq i \leq N,\] aunque en algunos casos resulta más conveniente utilizar \(r = sum_{i=1}^N r_i\). De hecho, veremos que ambos resultados son muy similares y podremos usar la suma en todos los casos.

Algoritmo para obtener \(P\):

  1. Fijar \(R\), \(t\), y \(0 < \epsilon < 1\). Por defecto fijamos \(\epsilon = 0.00001\).
  2. Obtener r.
  3. Calcular \(\hat{P}\).
  4. Calcular \(A = \hat{P}\); \(c = e^{rt}\); \(B = e^{rt}I\); \(sum = c\); \(k=1.\)
  5. Mientras que \(sum < 1-\epsilon\), calcular:

\[c = c*(rt)/k; \quad B = B + cA; \quad A = A \hat{P}\] \[sum = sum + c; \quad k = k + 1\] Al finalizar la matriz \(B\) es una aproximación de \(P(t)\) con error inferior a \(\epsilon\).

En realidad en nuestro algoritmo añadiremos un parámetro extra para indicar como se debe calcular el valor de \(r\), bien como el máximo los elementos por filas de \(R\) o como la suma de todos ellos.

matriz.prob.trans<- function(Rmat, ts, cal)
{
  # Algortimo de uniformización para obtener P(t)
  ################################################
  
  # Parámetros de la función
  # Rmat: matriz de tasas
  # ts: instante de tiempo
  # epsilon: error en la aproximación
  # cal: forms de obtener r con dos valores 1 = máximo, 2 = suma
  epsilon <- 1e-05
  # Paso 2. Calculo de r
  ris <- apply(Rmat, 1, sum)
  ifelse(cal == 1, rlimit <- max(ris), rlimit <- sum(Rmat))
  # Paso 3. Calculo de hat(P)
  hatP <- Rmat/rlimit
  diag(hatP) <- 1 - ris/rlimit
  # Paso 4. Calculo de matrices y vectores accesorios
  rts <- rlimit*ts
  A <- hatP
  c <- exp(-rts)
  B <- c*diag(nrow(Rmat))
  suma <- c
  k <- 1
  # Bucle simulación
  cota <- 1- epsilon
  while(suma < cota)
  {
    c <- c*rts/k
    B <- B + c*A
    A <- A%*%hatP
    suma <- suma + c
    k <- k + 1
  }
  return(round(B, 4))
}

Veamos un ejemplo donde la matriz de tasas viene dada por:

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

Si estamos interesados en las probabilidades de transición entre estados cuando han transcurrido 2 unidades de tiempo podemos calcular (ambas versiones):

Pmat1<-matriz.prob.trans(R, 2, 1); Pmat1
##        [,1]   [,2] [,3]   [,4]
## [1,] 0.2001 0.2001  0.4 0.1998
## [2,] 0.2002 0.2001  0.4 0.1997
## [3,] 0.1999 0.2000  0.4 0.2001
## [4,] 0.1998 0.1999  0.4 0.2003
Pmat2<-matriz.prob.trans(R, 2, 2); Pmat2
##        [,1]   [,2] [,3]   [,4]
## [1,] 0.2001 0.2001  0.4 0.1998
## [2,] 0.2002 0.2001  0.4 0.1997
## [3,] 0.1999 0.2000  0.4 0.2001
## [4,] 0.1998 0.1999  0.4 0.2003

Ambas matrices proporciona un resultado prácticamente idéntico. Estas matrices nos permiten obtener las probabildiades de pasar del estado 1 a cualquiera de los otros estados en dos unidades de tiempo, sin más que tomar los elementos de la fila 1 de la matriz obtenida.

Pmat1[1,]; Pmat2[1,]
## [1] 0.2001 0.2001 0.4000 0.1998
## [1] 0.2001 0.2001 0.4000 0.1998

¿Cómo interpretamos esas probabilidades?

Veamos ahora la aplicación de este algoritmo a alguno de los ejemplos con los que hemos ido trabajando en esta unidad.

Ejemplo 4.3 Para los datos correspondientes al cajero bancario estamos interesados en conocer la probabilidad de que después de 50 minutos de funcionamiento el sistema este completamente ocupado (1 usuario atendido y tres en cola), cuando partímos de \(0\) clientes en el sistema (\(p_{04}(4)\)). La matriz de tasas viene dada por:

estados <- c(0, 1, 2, 3, 4)
nestados <- length(estados)

R <- matrix(nrow = nestados, ncol = nestados, data = 0)
lambda <- 15/60
mu <- 1/6 

R[1,2] <- lambda 
R[2,1] <- mu 
R[2,3] <- lambda 
R[3,2] <- mu 
R[3,4] <- lambda 
R[4,3] <- mu 
R[4,5] <- lambda
R[5,4] <- mu

Obtenemos la distribución de probabilidad asociada al estado 4 para \(t = 50\):

# Matriz de probabilidades de transición
Pmat<-matriz.prob.trans(R, 50, 2)
Pmat
##        [,1]   [,2]   [,3]   [,4]   [,5]
## [1,] 0.0812 0.1190 0.1730 0.2528 0.3741
## [2,] 0.0793 0.1172 0.1722 0.2539 0.3775
## [3,] 0.0769 0.1148 0.1711 0.2553 0.3819
## [4,] 0.0749 0.1128 0.1702 0.2565 0.3856
## [5,] 0.0739 0.1118 0.1698 0.2571 0.3874

La probabilidad de interés es 0.3741 (\(p_{15}(50)\)) lo que demuestra que es factible que el sistema llegue al estado 4 partiendo del estado 0 después de 4 horas.

Aunque el cálculo teórico es muy preciso hay situaciones donde los sistemas reales con los que estamos trabajndo hacen bastante costoso obtener la matriz \(R\), y resulta más sencillo tratar de aproximar las probabilidades de transición mediante simulación. Para ello basta con replicar el sistema de simulación un número los suficientemente grande ya aproximar mediante Monte-Carlo. En este caso deberíamos obtner el estado dle sistema después de 50 minutos y aproximar la probabilidad como el número de veces que se alcanza el estado de interés dividido por las réplicas realizadas. la precisión de esta aproximación depende en gran medida del número de réplicas.

replicas <- 2500
envs <- lapply(1:replicas, function(i) {
    repcajero <- cola.MM1K(50, 15/60, 1/6, 1, 4)
})
# almacenamos análisis de recursos del sistema
simresource <- as_tibble(get_mon_resources(envs))
# Almacenamos el estado final de la cola en el último instante del sistema
salida <- simresource %>%
  group_by(replication) %>%
  summarise(estado = last(system))
# Estimamos la probabilidad
round(table(salida$estado)/replicas, 3)
## 
##     0     1     2     3     4 
## 0.090 0.123 0.163 0.250 0.374

Podemos ver que la aproximación obtenida es similar (hasta el segundo decimal) a la obtneida a partir de la matriz \(R\) del proceso. Simulando podemos aproximar las cantidades de interés siempre que el sistema empieza desde el punto de interés.

Para practicar la obtención de la matriz de probabilidades de transición puedes resolver los ejercicios B-1 a B-4 de la colección al final de la unidad.