Licensed under a Creative Commons BY-NC-SA 4.0 License.

1 Discrete-time Markov chains

1.1 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

1.2 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)

1.3 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)

1.4 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.

2 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.

2.1 Simmer: a simulation tool for R

Discrete Event Simulation is a method for simulating the behaviour and performance of real-life processes or systems.For example:

  • A supermarket has 4 check-outs. Customers arrive at the less crowded queue at rate \(\lambda\). Each customer is served within a mean time of \(1/\mu\).

Basically, using DES we will be able to model problems with the following characteristics:

  • The problem follow a continuous timeline which contains multiple discrete events.

  • All events can be classified as arrivals or departures from the system, having some processing in between.

This description defines a limited subset of all the problems that can be addressed following a Discrete-Event Simulation (DES) methodology but includes all those problems we want to address during this course.

In R, we have the package simmer. Besides the simulation and monitorization framework itself, the simmer.plot package offers convenience plotting functions based on the popular ggplot2 package.

library(simmer)
## Warning: package 'simmer' was built under R version 3.6.3
library(simmer.plot)
## Warning: package 'simmer.plot' was built under R version 3.6.3

2.1.1 Simmer: Trajectory

Simulating with simmer leverages the concept of trajectory. A trajectory is the path that a certain type of arrivals follows during their lifetime. Let’s take the previous supermarket example. The customer’s trajectory could be modelled as follows:

  1. Arrive at a check-out.
  2. Seize it (waiting in the queue if it is occupied).
  3. Wait an exponential service time \(\mu=4\).
  4. Release the check-out and leave.
# %>% operator concatenates operations

customer <- trajectory() %>%
  seize("check-out", amount=1) %>%
  timeout(function()rexp(1, 4)) %>%
  release("check-out", amount=1)

2.1.2 Simmer: Simulation

We can create a simulation environment with the necessary resources (4 check-outs, infinite queue) and a customer generator with rate \(\lambda=10\) following the previous trajectory:

supermarket <- simmer() %>%
  add_resource("check-out", capacity=4, queue_size=Inf) %>%
  add_generator("customer", customer, function() rexp(100, 10)) %>%
  run(until=20) # Simulation until the instant t=20

2.1.3 Simmer: Results visualization

We can visualise the average number of customers in the system over time using the following function:

# Store monitored data about arrivals, attributes and resources
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           1
# Plot results taking system column as item where the resource is called ckeck-out
## Adding steps = TRUE we will show the number of busy resources we have
plot(df.res, metric="usage", "check-out", items="system", steps=TRUE)

Horizontal dashed lines represent the limits for the queue, server and system (queue + server).We may add the option steps=TRUE to visualise the instantaneous absolute values.

2.2 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):

    1. 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)
    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)

3 Solving problems with R

3.1 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

3.2 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

3.3 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