A few weeks ago I wrote a tutorial on Markov Chains, where the example I gave was a model that generated text. I enjoyed writing that, but I realize if you’re interested in using Markov Chains for a more practical purpose, this tutorial only introduced the idea. So now, I want to provide an example that provides easier access to markov chain related programming.

Markov Chain – a random process, where we assume the previous state holds predictive power in predicting the next.

## Car Rental Example

Suppose a car rental agency has three locations in Ottawa: Downtown location, East end location and a West end location. The agency has a group of delivery drivers to serve all three locations. The agency’s statistician has determined the following:

- Of the calls to the Downtown location, 30% are delivered in Downtown area, 30% are delivered in the East end, and 40% are delivered in the West end.
- Of the calls to the East end location, 40% are delivered in Downtown area, 40% are delivered in the East end, and 20% are delivered in the West end
- Of the calls to the West end location, 50% are delivered in Downtown area, 30% are delivered in the East end, and 20% are delivered in the West end.

After making a delivery, a driver goes to the nearest location to make the next delivery. This way, the location of a specific driver is determined only by his or her previous location.

This is a markov chain because the current state has predictive power over the next state. The **state space **in this example is Downtown, East, and West.

Note – the code is also on my github.

We can model the markov chain by a transition matrix.

1 2 3 4 5 6 |
rentalStates <- c("Downtown", "East", "West") rentalTransition <- matrix(c(0.3, 0.3, 0.4, + 0.4, 0.4, 0.2, + 0.5, 0.3, 0.2), + byrow = T, nrow = 3, dimnames = list(rentalStates, rentalStates)) rentalTransition |

1 2 3 4 |
Downtown East West Downtown 0.3 0.3 0.4 East 0.4 0.4 0.2 West 0.5 0.3 0.2 |

There is a package in R that will make our lives easier. It’s called ‘markovchain’. Let’s define our model as a markovchain object. Calling the object is equivalent to it’s show() method.

We can even access elements of the transition matrix by indexing the object.

1 2 3 4 5 6 7 |
mcRental <- new("markovchain", states = rentalStates, byrow = T, transitionMatrix = rentalTransition, name = "Rental Cars") # We can access the transition matrix just by calling the mc object mcRental |

1 2 3 4 5 6 7 8 |
Rental Cars A 3 - dimensional discrete Markov Chain with following states Downtown East West The transition matrix (by rows) is defined as follows Downtown East West Downtown 0.3 0.3 0.4 East 0.4 0.4 0.2 West 0.5 0.3 0.2 |

You can also call a specific row of the matrix, using the mcRental object.

1 |
mcRental[1] # the probabilities that we go Downtown, East, and West, given that we are currently Downtown |

1 2 |
Downtown East West 0.3 0.3 0.4 |

We can plot it using plot(mcRental).

Ok, let’s ask some more statistical questions that might pop up if this were a real situation. Given that we are downtown, what is the probability being back downtown in exactly 2 trips?

To answer this, we have to find the each possible way we can arrive back downtown in 2 trips, and add their probabilities together.

1 2 3 4 5 6 7 8 9 10 11 12 |
# Here is a question to set up some of the functions # Given we are downtown, what is the probability we will be downtown in two trips? # We can go Downtown -> Downtown, a <- 0.3 * 0.3 # East -> Downtown (note that to we have to get the probability of going Downtown from the East location), b <- 0.3 * 0.4 # West -> Downtown (same logic here) c <- 0.4 * 0.5 a + b + c # The probability that we will be downtown in 2 trips. |

[1] 0.41

It turns out though, that you can get the same result by raising the transition matrix to the power n, where n is the number of trips.

1 |
mcRental ^ 2 |

1 2 3 4 5 6 7 8 |
Rental Cars^2 A 3 - dimensional discrete Markov Chain with following states Downtown East West The transition matrix (by rows) is defined as follows Downtown East West Downtown 0.41 0.33 0.26 East 0.38 0.34 0.28 West 0.37 0.33 0.30 |

The matrix shown above represents the probability you will go from one point to another in 2 trips.

We can do this for any number of trips. It’s makes logical sense that as the number of trips n increases, where you started will have less and less predictive power as where you end up. We can demonstrate this fact by raising the transition matrix to a sufficiently large number.

1 |
mcRental^20 |

1 2 3 4 5 6 7 8 |
Rental Cars^20 A 3 - dimensional discrete Markov Chain with following states Downtown East West The transition matrix (by rows) is defined as follows Downtown East West Downtown 0.3888889 0.3333333 0.2777778 East 0.3888889 0.3333333 0.2777778 West 0.3888889 0.3333333 0.2777778 |

This distribution, when the starting location is completely irrelevant, is known as the **stationary distribution**. It can be calculated with some linear algebra fanciness, but the markovchain package has a method called steadyStates(). Here is a question to make use of the method.

At the rental car company there are 70 drivers. Each driver is required to make 30 trips each day—how many drivers will end their day at the West location?

First we need to assure that 30 is sufficiently large, and we can do this by raising the transition matrix to 30 and assuring that there is no difference based on starting location.

After that, we can just multiply the scalar 70 times the stationary distribution, and that will give the vector of expected values.

1 2 |
mcRental ^ 30 70*steadyStates(mcRental) |

1 2 |
Downtown East West [1,] 27.22222 23.33333 19.44444 |

There are several more awesome functions in the markovchain package, two of which I want to highlight.

The conditional distribution just pulls out the row of the transition matrix. I love that there is a method called this, because it highlights the fact that markov chains have at their core conditional probability.

The summary function provides a more comprehensive view of the object. Each aspect of it (transient classes, irreducibility) are beyond the scope of this tutorial, but they’re worth noting.

1 |
summary(mcRental) |

1 2 3 4 5 6 7 8 9 |
Rental Cars Markov chain that is composed by: Closed classes: Downtown East West Recurrent classes: {Downtown,East,West} Transient classes: NONE The Markov chain is irreducible The absorbing states are: NONE |

1 |
conditionalDistribution(mcRental, "Downtown") |

1 2 |
Downtown East West 0.3 0.3 0.4 |

Sources –

- http://www.mast.queensu.ca/~stat455/lecturenotes/set3.pdf
- https://cran.r-project.org/web/packages/markovchain/vignettes/an_introduction_to_markovchain_package.pdf