Teoría de Redes - Práctica 3
Iñaki Úcar, Pablo Serrano, José Alberto Hernández, Ginés García Avilés
5 de marzo de 2018
Licenciado bajo Creative Commons BY-NC-SA 4.0.
1 Cadenas de Markov de tiempo discreto
1.1 Introducción
Para la manipulación de cadenas de Markov en R, vamos a utilizar la librería markovchain, creada por Giorgio Alfredo Spedicato. A continuación, nos familiarizaremos con su uso a través del siguiente modelo simple de predicción del tiempo visto en teoría (C: calor, F: frío):
\[\xymatrix{ *=<5mm,5mm>[o][F]{C} \ar@(ul,dl)[]_{0.6} \ar@/^/[r]^{0.4} & *=<5mm,5mm>[o][F]{F} \ar@(dr,ur)[]_{0.7} \ar@/^/[l]^{0.3} }\qquad\qquad P = \begin{pmatrix} 0.6 & 0.4 \\ 0.3 & 0.7 \end{pmatrix}\]
Esta introducción está basada en la documentación oficial.
En primer lugar, debemos importar la librería para tener acceso a sus clases y métodos. La siguiente instrucción equivale al #include de C o el import de Python:
library(markovchain)## Warning: package 'markovchain' was built under R version 3.6.31.2 Matriz de transición
Definimos la matriz de transición del ejemplo de arriba. Una matriz en R se define mediante el método matrix (véase la ayuda). Como entrada, recibe un vector de datos unidimensional, por lo que hay que especificar al menos una dimensión (número de filas o columnas) y cómo debe interpretar dicho vector (por filas o por columnas). En este caso, introducimos una matriz de 2x2 (ncol = 2) por filas (byrow = TRUE):
P <- matrix(data = c(0.6, 0.4,
0.3, 0.7),
ncol = 2, byrow = TRUE)1.3 Cadenas de Markov en R
Una cadena de Markov vendrá representada por un objeto de tipo “markovchain”. Para su creación, utilizamos el método genérico new e introducimos los nombres de los estados y la matriz de transición. Opcionalmente, podemos darle un nombre a la cadena, como se ve a continuación.
weather <- new("markovchain", name = "Simple weather model",
states = c("C", "F"), transitionMatrix = P)
weather## Simple weather model
## A 2 - dimensional discrete Markov Chain defined by the following states:
## C, F
## The transition matrix (by rows) is defined as follows:
## C F
## C 0.6 0.4
## F 0.3 0.7summary(weather)## Simple weather model Markov chain that is composed by:
## Closed classes:
## C F
## Recurrent classes:
## {C,F}
## Transient classes:
## NONE
## The Markov chain is irreducible
## The absorbing states are: NONEplot(weather)
En R no podemos utilizar los operadores habituales de multiplicación * y exponenciación ^ con matrices. La multiplicación se realiza con el siguiente operador:
P %*% P## [,1] [,2]
## [1,] 0.48 0.52
## [2,] 0.39 0.61Y para elevar una matriz a una potencia, necesitamos cargar la librería expm:
library(expm)## Warning: package 'expm' was built under R version 3.6.3## Loading required package: Matrix##
## Attaching package: 'expm'## The following object is masked from 'package:Matrix':
##
## expmP%^%4## [,1] [,2]
## [1,] 0.4332 0.5668
## [2,] 0.4251 0.5749En cambio, el objeto markovchain se comporta aritméticamente como una matriz de transición sobre la que sí podemos aplicar los operadores habituales. La siguiente tabla detalla las principales operaciones definidas sobre él.
| Método | Resultado |
|---|---|
mcA * mcB |
Multiplicación de matrices de transición |
mcA^4 |
Exponenciación de la matriz de transición |
mcA[1,] |
Acceso a la fila 1 de la matriz de transición |
mcA[, 3] |
Acceso a la columna 3 de la matriz de transición |
mcA[3, 2] |
Acceso al elemento (3, 2) de la matriz de transición |
t(mcA) |
Trasposición de la matriz de transición |
mcA == mcB |
Comparación entre dos objetos (estados y matriz) |
dim(mcA) |
Dimensión de la matriz de transición |
states(mcA) |
Extracción de los estados |
Por tanto, podemos aplicar este objeto directamente a la ecuación de Chapman-Kolmogorov:
\[\pi^{(n)} = \pi^{(0)}P^n\]
Partimos del estado inicial “F” y obtenemos el vector de distribución de probabilidad para \(n = 6\):
initialState <- c(0, 1)
newState <- initialState * weather^6
newState## C F
## [1,] 0.428259 0.571741Pero, además, el objeto markovchain nos permite obtener mucha información sobre la cadena: irreducibilidad, accesibilidad, clases, tipos de estados, etc.
is.irreducible(weather)## [1] TRUEis.accessible(weather, from = "C", to = "F")## [1] TRUEperiod(weather)## [1] 1absorbingStates(weather)## character(0)transientStates(weather)## character(0)steadyStates(weather)## C F
## [1,] 0.4285714 0.5714286Lo más interesante de esta librería es que nos permite realizar de forma sencilla simulación y estimación de cadenas de Markov. La simulación, dada una cadena de Markov previamente definida, se realiza de manera análoga a la simulación de variables aleatorias convencionales (rnorm, rexp, etc.) mediante la función rmarkovchain.
Generación de 1000 muestras aleatorias a partir de la cadena de Markov que hemos definido al comienzo del tutorial. El estado inicial es aleatorio:
x <- rmarkovchain(1000, weather)Generación de 1000 muestras aleatorias partiendo del estado inicial “C”:
y <- rmarkovchain(1000, weather, t0 = "C")La estimación consiste en ajustar unas probabilidades de transición a partir de un vector de muestras aleatorias. El ajuste a partir del vector x generado anteriormente se aproxima bastante a la cadena original:
weatherFit <- markovchainFit(x)
plot(weatherFit$estimate)
1.4 Convergencia
Dada una cadena de Markov irreducible y aperiódica, se puede demostrar que existe un único estado estacionario al que la cadena converge cuando \(n\to\infty\) independientemente del estado inicial. Tomando nuestro modelo de predicción del tiempo, podemos ilustrar esta propiedad elevando la matriz de transición a una potencia suficientemente alta y observando que todas las filas de la matriz resultante tienden al estado estacionario:
is.irreducible(weather)## [1] TRUEweather^100## Simple weather model^100
## A 2 - dimensional discrete Markov Chain defined by the following states:
## C, F
## The transition matrix (by rows) is defined as follows:
## C F
## C 0.4285714 0.5714286
## F 0.4285714 0.5714286steadyStates(weather)## C F
## [1,] 0.4285714 0.5714286Esto es equivalente a decir que, cuando simulamos una cadena de Markov (obtenemos un conjunto de muestras), las probabilidades de los estados serán estables al alcanzar un número de muestras lo suficientemente elevado, y no van a cambiar por añadir más.
samples <- rmarkovchain(1000, weather)
# Proporción de cada estado acumulada en cada paso
yC <- cumsum(samples=="C") / seq_along(samples)
yF <- cumsum(samples=="F") / seq_along(samples)
# Evolución de la distribución de C
plot(yC, type="l", col="red", ylim=c(0,1))
# Evolución de la distribución de F
lines(yF, col="green")
# Valores estacionarios teóricos
abline(steadyStates(weather)[1], 0, lty=2, col="red")
abline(steadyStates(weather)[2], 0, lty=2, col="green")
En la figura anterior se aprecia que la simulación converge. Pero, ¿qué potencia es lo suficientemente alta, qué número de muestras el lo suficientemente elevado para asegurar una buena aproximación al estado estacionario? Aquí entra en juego el llamado tiempo de convergencia, que nos dará una idea de a qué velocidad se va a reducir el error conforme avanza el tiempo.
Podemos definir el error, de forma simple, como el sumatorio de las diferencias absolutas entre componentes de dos potencias sucesivas de la matriz de transición, es decir,
\[\epsilon(n) = \sum_i\sum_j|\left(P^n\right)_{ij} - \left(P^{n-1}\right)_{ij}|\]
Vamos a construir una función en R para obtener este error y poder visualizarlo.
# Función de error
# @param n es el tiempo (valor o conjunto de valores)
# @param mc es la cadena de Markov
err <- function(n, mc) {
# Comprobamos que los valores de entrada son números naturales
if (!isTRUE(all(n > 0 && n == floor(n))))
stop("'n' must only contain positive integer values")
# Reservamos memoria para la salida
res <- numeric(length(n))
# Para cada valor de entrada "n", calculamos err(n)
for (i in 1:length(n)) {
Pn <- (mc^n[i])@transitionMatrix
Pn1 <- (mc^(n[i]-1))@transitionMatrix
res[i] <- sum(abs(Pn - Pn1))
}
return(res)
}
# Representación de la función de error
x <- 1:10
y <- err(x, weather)
plot(x, y, type="o")
Comprobamos que el error decae exponencialmente.
2 Continuous-time Markov chains
Ejemplo 1: Sea una gasolinera con un solo surtidor y sin espacio para que los vehículos esperen (si un vehículo llega a la gasolinera y el surtidor está ocupado, se va). A la gasolinera llegan vehículos según un proceso de Poisson de tasa \(\lambda\) igual a 3 cada 20 minutos, de los cuales \(prob(c)=\) 75% son coches y \(prob(m)=\) 25% motocicletas. El tiempo que le lleva a un usuario dejar libre el surtidor (es decir, repostar y pagar) se puede modelar como una variable aleatoria exponencial, de media \(1/\mu_c = 8\) minutos si es un coche, y \(1/\mu_m=3\) minutos si es una moto.
Este problema se puede modelar con la siguiente cadena de Markov de tiempo continuo:
\[\xymatrix{ *=<15mm,8mm>[o][F]{car} \ar@/^/[r]^{\mu_c} & *=<15mm,8mm>[o][F]{empty} \ar@/^/[l]^{\lambda \cdot prob(c)} \ar@/^/[r]^{\lambda \cdot prob(m)} & *=<15mm,8mm>[o][F]{m/cycle} \ar@/^/[l]^{\mu_m} }\qquad\qquad Q = \begin{pmatrix} -1/8 & 1/8 & 0 \\ 0.75\cdot 3/20 & -3/20 & 0.25\cdot 3/20 \\ 0 & 1/3 & -1/3 \end{pmatrix}\]
La cadena es irreducible y recurrente. Para encontrar la distribución de estado estacionaria teóricamente, debemos resolver las ecuaciones de balance
\[pQ = 0\]
a las que hay que añadir la relación de cierre
\[\sum_i p_i = 1\]
para obtener un sistema compatible determinado. Dado que hay \(\operatorname{dim}(Q)-1\) columnas independientes, dicha relación de cierre se traduce en sustituir cualquier columna por unos e igualar a uno dicha componente al otro lado de la ecuación, es decir:
\[p\begin{pmatrix} 1 & 1/8 & 0 \\ 1 & -3/20 & 0.25\cdot 3/20 \\ 1 & 1/3 & -1/3 \end{pmatrix} = (1, 0, 0)\]
Al resolver el sistema, obtendremos un vector \(p\) que representa la probabilidad de estar en cada uno de los estados a largo plazo. Por tanto, para calcular el número medio de usuarios, estamos interesados en sumar las probabilidades multiplicadas por el número de usuarios en sus respectivos estados. En este caso, \(1\cdot p_1 + 0\cdot p_2 + 1\cdot p_3\).
lambda <- 3/20 # Tasa de llegadas
mu <- c(1/8, 1/3) # Tasa de servicio de coches y motos respectivamente
p <- 0.75 # Probabilidad de que el vehículo sea coche
A <- matrix(c(1, mu[1], 0,
1, -lambda, (1-p)*lambda,
1, mu[2], -mu[2]), byrow=T, ncol=3)
P <- solve(t(A), c(1, 0, 0))
N_average_theor <- sum(P * c(1, 0, 1))
N_average_theor## [1] 0.5031056En el código anterior, la función t traspone la matriz para resolver el sistema con solve. Consulte la ayuda para más detalles.
2.1 Simmer: una herramienta de simulación para R
La simulación de eventos discretos (DES) es un método para simular comportamiento y rendimiento de procesos/sistemas del mundo real. Por ejemplo:
Supongamos que queremos simular por ordenador la evolución de un problema como los siguientes:
- Un supermercado tiene en todo momento 4 cajas en funcionamiento. Los clientes se van colocando en la cola menos ocupada según van llegando a una tasa \(\lambda\). Cada cliente es atendido en un tiempo medio \(1/\mu\).
Básicamente, usando DES seremos capaces de modelar problemas con las siguientes características:
Siguen una línea de tiempo continua en la que se produce una secuencia de eventos discretos de distinto tipo.
Todos los eventos se pueden clasificar en llegadas o salidas del sistema, entre las cuales se produce un procesamiento.
Esta descripción responde a un subconjunto limitado, pero de interés para esta asignatura, de todos los problemas que pueden abordarse mediante una metodología de simulación por eventos discretos (DES).
En R, disponemos del paquete simmer. Además del framework de simulación y monitorización de eventos, el paquete simmer.plot nos ofrece cómodas funciones de representación basadas en el popular paquete ggplot2.
library(simmer)## Warning: package 'simmer' was built under R version 3.6.3library(simmer.plot)## Warning: package 'simmer.plot' was built under R version 3.6.32.1.1 Simmer: Trayectoria
La simulación con simmer se basa en el concepto de trayectoria. Una trayectoria es el camino que cierta clase de llegadas van a seguir durante su tiempo de vida. Tomemos el ejemplo anterior del supermercado. La trayectoria para los clientes podría modelarse como sigue:
- Llegada a la caja.
- Acceso a la misma (con espera en cola si se encuentra ocupada).
- Espera de un tiempo exponencial correspondiente al tiempo de servicio \(\mu=4\).
- Liberación de la caja.
Y en R (nótese la conveniencia del operador pipe %>%, que se encarga de concatenar acciones una detrás de otra):
customer <- trajectory() %>%
seize("check-out", amount=1) %>%
timeout(function() rexp(1, 4)) %>%
release("check-out", amount=1)2.1.2 Simmer: Simulación
Podemos crear un entorno de simulación con el recurso necesario (4 cajas, con cola de espera infinita) y un generador de clientes con tasa \(\lambda=10\) que sigan la trayectoria anterior:
supermarket <- simmer() %>%
add_resource("check-out", capacity=4, queue_size=Inf) %>%
add_generator("customer", customer, function() rexp(100, 10)) %>%
run(until=20) # Simulación hasta el instante t=202.1.3 Simmer: Visualización de resultados
Podemos visualizar la ocupación media del sistema a lo largo del tiempo mediante la siguiente función:
df.res <- get_mon_resources(supermarket)
head(df.res)## resource time server queue capacity queue_size system limit
## 1 check-out 0.04516067 1 0 4 Inf 1 Inf
## 2 check-out 0.08772235 2 0 4 Inf 2 Inf
## 3 check-out 0.14039323 3 0 4 Inf 3 Inf
## 4 check-out 0.14761076 2 0 4 Inf 2 Inf
## 5 check-out 0.15066620 1 0 4 Inf 1 Inf
## 6 check-out 0.18132978 2 0 4 Inf 2 Inf
## replication
## 1 1
## 2 1
## 3 1
## 4 1
## 5 1
## 6 1plot(df.res, metric="usage", "check-out", items="system", steps=TRUE)
Las líneas horizontales discontinuas representan los límites de la cola, el servidor y el sistema (cola + servidor).Si además estamos interesados en visualizar la ocupación instantánea, añadimos la opción steps=TRUE
2.2 Convergencia
Siguiendo con el Ejemplo 1, tras obtener la solución teórica, procedemos a simular el sistema y comprobar que converge a la solución teórica. Para ello, usaremos simulación de eventos discretos (DES):
- Dividimos el problema en dos trayectorias, porque los tiempos de servicio son distintos para el coche y la moto.
car <- trajectory() %>%
seize("pump", amount=1) %>%
timeout(function() rexp(1, mu[1])) %>%
release("pump", amount=1)
mcycle <- trajectory() %>%
seize("pump", amount=1) %>%
timeout(function() rexp(1, mu[2])) %>%
release("pump", amount=1)- Definimos un simulador que contiene dos generadores, uno por cada tipo de vehículo que llega y el recurso compartido que tiene la gasolinera (distribuidor de gasolina).
gas.station <- simmer() %>%
# Un solo recurso sin cola
add_resource("pump", capacity=1, queue_size=0) %>%
# Generador para los coches de tasa = p * lambda
add_generator("car", car, function() rexp(100, p*lambda)) %>%
# Generador para las motos de tasa = (1 - p) * lambda
add_generator("mcycle", mcycle, function() rexp(100, (1-p)*lambda)) %>%
run(1000/lambda)
# Evolución del número medio de usuarios en el sistema + Valor teórico
plot(get_mon_resources(gas.station), metric="usage", "pump", items="system") +
geom_hline(yintercept = N_average_theor)
3 Resolución de problemas con R
3.1 Problema 1
El problema del gato y el ratón. Hay cinco cajas en fila, con un gato en la primera y un ratón en la última. En cada intervalo de tiempo \(T\), cada animal salta a una caja colindante, completamente al azar. El juego termina cuando ambos animales saltan a la misma caja. Calcule la duración media del juego.
P <- matrix(data = c(0, 1, 0, 0, 0,
0.25, 0, 0.25, 0.25, 0.25,
0, 0.5, 0, 0, 0.5,
0, 0.5, 0, 0, 0.5,
0, 0, 0, 0, 1),
ncol = 5, byrow = TRUE)
catmouse <- new("markovchain", name = "Cat & Mouse problem",
states = c("1,5", "2,4", "1,3", "3,5", "End"), transitionMatrix = P)
n15End <- replicate(1000, match("End", rmarkovchain(100, catmouse, t0 = "1,5")))
mean(n15End) # ~4.5## [1] 4.3973.2 Problema 2
Una tarjeta de red tiene tres velocidades de transmisión, 1, 4 y 8 Mbps, que dan lugar a una probabilidad de pérdida de trama de \(p_1=1/2\), \(p_4=1/2\) y \(p_8=1/4\) respectivamente. El esquema que se usa para adaptar la velocidad de transmisión es el siguiente: ante una pérdida, se transmite a 1 Mbps, mientras que tras dos éxitos consecutivos a una velocidad se aumenta la velocidad al siguiente valor superior disponible. Calcule la velocidad media (en Mbps) de transmisión.
#{r, echo=FALSE, eval=FALSE}
P <- matrix(
c(1/2, 1/2, 0, 0, 0,
1/2, 0, 1/2, 0, 0,
1/2, 0, 0, 1/2, 0,
1/2, 0, 0, 0, 1/2,
1/4, 0, 0, 0, 3/4),
ncol = 5, byrow = TRUE)
chain <- new("markovchain",
states = c("1.1", "1.2", "4.1", "4.2", "8"),
transitionMatrix = P)
plot(chain)
sum(steadyStates(chain) * c(1, 1, 4, 4, 8))## [1] 2.3529413.3 Problema 3
Sea un edificio con dos ascensores, cada uno con un tiempo de vida que se puede modelar con una variable aleatoria exponencial de media \(1/\lambda\). Hay dos políticas de reparación:
- Se avisa a mantenimiento cuando los dos ascensores están estropeados, siendo el tiempo de reparación de los dos ascensores (a la vez) una variable aleatoria exponencial de media \(1/\mu\).
- Se avisa a mantenimiento en cuanto un ascensor se estropea, siendo el tiempo de reparación por ascensor una variable aleatoria exponencial de media \(2/\mu\), y sólo se repara un ascensor a cada vez.
Calcule la proporción de tiempo que los dos ascensores están funcionando, si \(1/\lambda = 2\)~semanas y \(1/\mu= 1\)~semana.
#{r, echo=FALSE, eval=FALSE}
lambda <- 1/2
mu <- 1
### a) ###
##########
Q <- matrix(
c(1, 2*lambda, 0,
1, -lambda, lambda,
1, 0, -mu),
ncol = 3, byrow = TRUE)
# Estado inicial
P_a <- solve(t(Q), c(1, 0, 0))
# Seleccionamos el valor donde los dos ascensores estan activos
P_a[1]## [1] 0.25### b) ###
##########
Q <- matrix(
c(1, 2*lambda, 0,
1, -lambda-(mu/2), lambda,
1, mu/2, -mu/2),
ncol = 3, byrow = TRUE)
# Estado inicial
P_b <- solve(t(Q), c(1, 0, 0))
# Seleccionamos el valor donde los dos ascensores estan activos
P_b[1]## [1] 0.2A continuación, calculamos los resultados usando simulaciones:
lambda <- 1/2
mu <- 1
### a) ###
##########
FAIL <- function() rexp(1, lambda)
REPAIR <- function() rexp(1, mu)
elevator <- trajectory() %>%
seize("working") %>%
timeout(FAIL) %>%
release("working") %>%
batch(2) %>%
timeout(REPAIR) %>%
separate() %>%
rollback(6)
building <- simmer() %>%
add_resource("working", capacity=Inf) %>%
add_generator("elevator", elevator, at(0, 0)) %>%
run(40000/lambda)
building.df.res <- get_mon_resources(building)
# Obtener el array de marcas de tiempo
deltas <- diff(building.df.res$time)
# Calculamos las posiciones del array de tiempos donde los dos elevadores están trabajando
both_working <- which(building.df.res$system == 2)
t_both_working <- deltas[both_working]
# Calcualmos el tiempo dividiendo el tiempo total en el que los dos ascensores estan trabajando entre el tiempo total
sum(t_both_working, na.rm=TRUE) / max(building.df.res$time)## [1] 0.2501716P_a[1]## [1] 0.25### b) ###
##########
FAIL <- function() rexp(1, lambda)
REPAIR <- function() rexp(1, mu/2)
elevator_2 <- trajectory() %>%
seize("working_2") %>%
timeout(FAIL) %>%
release("working_2") %>%
seize("repairperson") %>%
timeout(REPAIR) %>%
release("repairperson") %>%
rollback(6)
building_2 <- simmer() %>%
add_resource("working_2", capacity=Inf) %>%
add_resource("repairperson", capacity=1, queue_size=Inf) %>%
add_generator("elevator_2", elevator_2, at(0, 0)) %>%
run(40000/lambda)
# Como tenemos dos recursos (working_2 y repariperson), seleccionamos los resultados del recurso en cuestion
building_2.df.res <- get_mon_resources(building_2)%>%
subset(resource == "working_2")
deltas_2 <- diff(building_2.df.res$time)
both_working_2 <- which(building_2.df.res$system == 2)
t_both_working_2 <- deltas_2[both_working_2]
sum(t_both_working_2, na.rm=TRUE) / max(building_2.df.res$time)## [1] 0.202305P_b[1]## [1] 0.2