Python Until Time T 100 Simulate a Continuous time Markov Chain With Generator
Licensed under a Creative Commons BY-NC-SA 4.0 License.
Discrete-time Markov chains
Introduction
We plan to use the markovchain library, by Giorgio Alfredo Spedicato, to deal with Markov chains in R. This first section is devoted to introduce its classes and methods using a simple weather prediction model (W: warm, C: cold):
\[\xymatrix{ *=<5mm,5mm>[o][F]{W} \ar@(ul,dl)[]_{0.6} \ar@/^/[r]^{0.4} & *=<5mm,5mm>[o][F]{C} \ar@(dr,ur)[]_{0.7} \ar@/^/[l]^{0.3} }\qquad\qquad P = \begin{pmatrix} 0.6 & 0.4 \\ 0.3 & 0.7 \end{pmatrix}\]
This intro was based on the official documentation.
First of all, we need to import the new library to get access to it. The following line is equivalent to C's #include or Python's import:
library(markovchain) ## Warning: package 'markovchain' was built under R version 3.6.3 Transition matrix definition
To start with, we define the transition matrix of the example above using the function matrix (see ?matrix). The components are arranged in a plain vector, so we must specify one dimension at least (number of rows or columns) and tell how to read this vector (by rows or columns). In our example, we define a 2x2 matrix (ncol = 2) by rows (byrow = TRUE):
P <- matrix(data = c(0.6, 0.4, 0.3, 0.7), ncol = 2, byrow = TRUE) Markov chains in R
A Markov chain is described by a "markovchain" object. We'll use the function new to create this type of objects by feeding it with the states' name and the transition matrix. Optionally, we can set a name for the chain.
weather <- new("markovchain", name = "Simple weather model", states = c("W", "C"), transitionMatrix = P) weather ## Simple weather model ## A 2 - dimensional discrete Markov Chain defined by the following states: ## W, C ## The transition matrix (by rows) is defined as follows: ## W C ## W 0.6 0.4 ## C 0.3 0.7 summary(weather) ## Simple weather model Markov chain that is composed by: ## Closed classes: ## W C ## Recurrent classes: ## {W,C} ## Transient classes: ## NONE ## The Markov chain is irreducible ## The absorbing states are: NONE plot(weather) 
We cannot use the common multiplication * and exponentiation ^ operators with R matrices. The product of two matrices can be done as follows:
P %*% P ## [,1] [,2] ## [1,] 0.48 0.52 ## [2,] 0.39 0.61 And the matrix exponentiation is defined in the non-standard library 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': ## ## expm P%^%4 ## [,1] [,2] ## [1,] 0.4332 0.5668 ## [2,] 0.4251 0.5749 However, a markovchain behaves like a transition matrix and accepts the common arithmetic operators. The following table lists the main operations defined over this object.
| Method | Result |
|---|---|
mcA * mcB | Transition matrix multiplication |
mcA^4 | Transition matrix exponentiation |
mcA[1,] | Accessing row 1 of the ransition matrix |
mcA[, 3] | Accessing column 3 of the ransition matrix |
mcA[3, 2] | Accessing element (3, 2) of the ransition matrix |
t(mcA) | Transition matrix transposition |
mcA == mcB | Comparing two objects (states and matrix) |
dim(mcA) | Transition matrix dimension |
states(mcA) | Extracting the states |
Therefore, we can use this object directly to operate with the Chapman-Kolmogorov equation:
\[\pi^{(n)} = \pi^{(0)}P^n\]
We start from the initial state "C" and obtain the probability distribution vector for \(n = 6\):
initialState <- c(0, 1) newState <- initialState * weather^6 newState ## W C ## [1,] 0.428259 0.571741 Moreover, the markovchain object gives us a lot of useful information about the chain: irreducibility, accessibility, classes, types of states and so on and so forth.
is.irreducible(weather) ## [1] TRUE is.accessible(weather, from = "W", to = "C") ## [1] TRUE period(weather) ## [1] 1 absorbingStates(weather) ## character(0) transientStates(weather) ## character(0) steadyStates(weather) ## W C ## [1,] 0.4285714 0.5714286 The more interesting part lies on the simulation and estimation capabilities of this library. Given a Markov chain, simulation is performed in the same way as conventional random variables (rnorm, rexp, etc.) using the function rmarkovchain.
Generation of 1000 random samples from the "weather" chain with random initial state:
x <- rmarkovchain(1000, weather) Generation of 1000 random samples with initial state "W":
y <- rmarkovchain(1000, weather, t0 = "W") Estimation consists of fitting the transition probabilities from a vector of random samples. Fitting from the previous x vector is close to the original chain:
weatherFit <- markovchainFit(x) plot(weatherFit$estimate) 
Convergence
Given an irreducible and aperiodic Markov chain, it can be shown that there exists a single steady state to which the chain converges when \(n\to\infty\) with independence of the initial state. Using our weather prediction model, let's illustrate this property by raising the transition matrix to a high enough power and observing that all rows of the resulting matrix tend to the steady state:
is.irreducible(weather) ## [1] TRUE weather^100 ## Simple weather model^100 ## A 2 - dimensional discrete Markov Chain defined by the following states: ## W, C ## The transition matrix (by rows) is defined as follows: ## W C ## W 0.4285714 0.5714286 ## C 0.4285714 0.5714286 steadyStates(weather) ## W C ## [1,] 0.4285714 0.5714286 What is that but to say, when simulating a Markov chain (obtaining random samples), that the states' ratio will be steady with a high enough number of samples.
samples <- rmarkovchain(1000, weather) # Cumulative ratios yW <- cumsum(samples=="W") / seq_along(samples) yC <- cumsum(samples=="C") / seq_along(samples) # Evolution of the distribution of "W" plot(yW, type="l", col="red", ylim=c(0,1)) # Evolution of the distribution of "C" lines(yC, col="green") # Theoretical values abline(steadyStates(weather)[1], 0, lty=2, col="red") abline(steadyStates(weather)[2], 0, lty=2, col="green") 
The above figure shows that the simulation converge. But, which power is high enough?, how many samples are sufficient to ensure a good approximation to the steady state? The key point here is the convergence time, which tells us the speed of error decay.
The error can be defined as the sum of the absolute differences between elements of two consecutive powers of the transition matrix, that is,
\[\epsilon(n) = \sum_i\sum_j|\left(P^n\right)_{ij} - \left(P^{n-1}\right)_{ij}|\]
Let's implement an R function to obtain this error.
# Error function # @param n is the time (one value or a vector of values) # @param mc is the Markov chain err <- function(n, mc) { # Check that the input values are positive integers if (!isTRUE(all(n > 0 && n == floor(n)))) stop("'n' must only contain positive integer values") # Reserve some memory for the result res <- numeric(length(n)) # For each n, calculate 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) } # Create an array of values x <- 1:10 # Call the previously created function to calculate the error y <- err(x, weather) plot(x, y, type="o") 
As can be seen, the error decays exponentially.
Continuous-time Markov Chains
Example 1: A gas station has a single pump and no space for vehicles to wait (if a vehicle arrives and the pump is not available, it leaves). Vehicles arrive to the gas station following a Poisson process with a rate \(\lambda\) of 3 every 20 minutes, of which \(prob(c)=\) 75% are cars and \(prob(m)=\) 25% are motorcycles. The refuelling time can be modelled with an exponential random variable with mean \(1/\mu_c = 8\) minutes for cars and \(1/\mu_m = 3\) minutes for motorcycles.
This problem is described by the following continuous-time Markov chain:
\[\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}\]
The chain is irreducible and recurrent. To theoretically find the steady state distribution, we have to solve the balance equations
\[pQ = 0\]
with the constraint
\[\sum_i p_i = 1\]
There are \(\operatorname{dim}(Q)-1\) independent columns, so the latter constraint is equivalent to substitute any column by ones and match it to one at the other side of the equation, that is:
\[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)\] The solution \(p\) represents the probability of being at each state in the long-term. Therefore, we can calculate the average number of customers in the system by summing these probabilities multiplied by the number of customers at each state. In our case, \(1\cdot p_1 + 0\cdot p_2 + 1\cdot p_3\).
lambda <- 3/20 # Arrival rate mu <- c(1/8, 1/3) # Service rate (cars, motorcycles) p <- 0.75 # Probability of car 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.5031056 In the previous chunk, the t function trasposes the matrix to solve the system of equations using solve. See the help for more details.
Convergence
Following with the Example 1, after obtaining the theoretical solution of the problem, we want to verify that the system exactly converges at that point. In order to do so, we are going to use discrete event simulation (DES):
-
- Break down the problem into two trajectories, because the service times are different.
# Car trajectory car <- trajectory() %>% seize("pump", amount=1) %>% timeout(function() rexp(1, mu[1])) %>% release("pump", amount=1) # Motorcycle trajectory mcycle <- trajectory() %>% seize("pump", amount=1) %>% timeout(function() rexp(1, mu[2])) %>% release("pump", amount=1) -
- Define a simulator with two generators, one for each type of vehicle and the shared resource that the gas station has (pump).
gas.station <- simmer() %>% # One server, no queue add_resource("pump", capacity=1, queue_size=0) %>% # Car generator of rate = p * lambda add_generator("car", car, function() rexp(100, p*lambda)) %>% # Motorcycle generator of rate = (1 - p) * lambda add_generator("mcycle", mcycle, function() rexp(100, (1-p)*lambda)) %>% run(1000/lambda) # Evolution of the average number of customers in the system + Theoretical value using ggplot horizontal line plot(get_mon_resources(gas.station), metric="usage", "pump", items="system") + geom_hline(yintercept = N_average_theor) 
Solving problems with R
Problem 1
The cat and mouse problem. Assume there are five boxes in a row, and there's a cat on the first one and a mouse on the last one. Each time interval \(T\), each animal jumps to a box nearby, which is chosen at random (unless they are on the first or last box). The "game" ends when they both jump into the same box, for obvious reasons. Compute the average length of the game.
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) P <- matrix(data = c(0, 1, 0, 0, 0.25, 0, 0.5, 0.25, 0, 0.5, 0, 0.5, 0, 0, 0, 1), ncol = 4, byrow = TRUE) catmouse <- new("markovchain", name = "Cat & Mouse problem", states = c("1,5", "2,4", "1,3,5", "End"), transitionMatrix = P) n15End <- replicate(100000, match("End", rmarkovchain(100, catmouse, t0 = "1,5"))) mean(n15End) # ~4.5 ## [1] 4.49942 Problem 2
A wireless interface has three transmission rates, 1, 4 and 8 Mbps, that result in three frame loss probabilities, \(p_1=1/2\), \(p_4=1/2\) and \(p_8=1/4\) respectively. The rate adaptation rate works as follows: whenever there is a loss, the transmission rate is set to 1 Mbps, while whenever there are two successess in a row, the transmission rate increases to the next available rate. Compute the average transmission rate (in Mbps).
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.352941 Problem 3
Assume a building with two elevators, each with a lifetime that can be modeled with an exponential random variable of average \(1/\lambda\). There are two repair policies:
a)Maintenance is notified when the two elevators are broken, with the repair time of the two elevators (at the same time) being an exponential random variable of average \(1/\mu\).
b)Maintenance is notified as soon as any elevator breaks, with the elevator repair time being an exponential random variable of average \(2/\mu\) (i.e., twice the former), and only one elevator is repaired at a time.
Calculate the proportion of time that the two elevators are running, if \(1/\lambda = 2\)~weeks and \(1/\mu = 1\)~week.
#{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) # Initial state 1, 0, 0 P_a <- solve(t(Q), c(1, 0, 0)) # Where both elevators are running 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) # Initial state 1, 0, 0 P_b <- solve(t(Q), c(1, 0, 0)) # Where both elevators are running P_b[1] ## [1] 0.2 We next compute this result using simulations:
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) # Retrieve an array of times deltas <- diff(building.df.res$time) # Compute array positions where both elevators were working both_working <- which(building.df.res$system == 2) t_both_working <- deltas[both_working] # Compute the time by dividing the time where we had 2 elevators working over the total time sum(t_both_working, na.rm=TRUE) / max(building.df.res$time) ## [1] 0.2465291 P_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) # Since we have now two resources, it is required to select only those related to working_2 building_2.df.res <- get_mon_resources(building_2)%>% subset(resource == "working_2") # Retrieve an array of timestamps deltas_2 <- diff(building_2.df.res$time) # Compute array positions where both elevators were working both_working_2 <- which(building_2.df.res$system == 2) t_both_working_2 <- deltas_2[both_working_2] # Compute the time by dividing the time where we had 2 elevators working over the total time sum(t_both_working_2, na.rm=TRUE) / max(building_2.df.res$time) ## [1] 0.2004387 P_b[1] ## [1] 0.2 Source: https://www.it.uc3m.es/gigarcia/index_files/labs/03-R-en.html
0 Response to "Python Until Time T 100 Simulate a Continuous time Markov Chain With Generator"
Post a Comment