---
title: "Markov Chains"
author: "Prof. Wright"
date: "April 16, 2018"
output: pdf_document
header-includes:
- \usepackage{tikz}
- \usetikzlibrary{arrows,automata}
- \usepackage{amsmath}
---
```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```
## Markov Chain Review
A Markov chain on $n$ states $S_1, S_2,\dots , S_n$ is defined by a $n\times n$ transition matrix $P$, were the entry $P_{i,j}$ in row $i$, column $j$ gives the probability $P[j\rightarrow i]$ of transitioning from state $S_j \rightarrow S_j$.
Here is a *state diagram* for three states.
\begin{center}
\begin{tikzpicture}[->,>=stealth,shorten >=1pt,auto,node distance=120pt,thick]
\tikzstyle{every state}=[fill=blue,draw=none,text=white,minimum size=30pt]
\tikzset{every edge/.append style={font=\small}}
\node[state] (A) {$S_1$};
\node[state] (B) [below left of=A] {$S_2$};
\node[state] (C) [below right of=A] {$S_3$};
\path
(A) edge [out=45,in=135,distance=40] node [pos=.5,above] {$P[1\rightarrow 1]$} (A)
(A) edge [bend left=20] node [pos=.4,left] {$P[1\rightarrow 2]$} (B)
(A) edge [bend left=30] node {$P[1\rightarrow 3]$} (C)
(B) edge [out=45+120,in=135+120,distance=40] node [pos=.5,left] {$P[2\rightarrow 2]$} (B)
(B) edge [bend left=30] node {$P[2\rightarrow 1]$} (A)
(B) edge [bend left=20] node [below] {$P[2\rightarrow 3]$} (C)
(C) edge [out=45-120,in=135-120,distance=40] node [pos=.5,right] {$P[3\rightarrow 3]$} (C)
(C) edge [bend left=20] node [pos=.6,right] {$P[3\rightarrow 1]$} (A)
(C) edge [bend left=20] node [pos=.5] {$P[3\rightarrow 2]$} (B);
\end{tikzpicture}
\end{center}
The transition matrix of probabilities, with $p_{i,j} = P[j\rightarrow i]$, is
$$ P = \begin{bmatrix}
p_{11} & p_{12}& p_{13}\\
p_{21} & p_{22}& p_{23}\\
p_{31} & p_{32}& p_{33}
\end{bmatrix} $$
Under very general hypotheses, the transition matrix $P$, will have a steady state vector $v_{ss}$ satisfying
$$Pv_{ss}=v_{ss}.$$
In other words, $v_{ss}$ is an eigenvector of $P$ corresponding to eigenvalue 1.
For example, define the following transition matrix of probabilities (note that each column sums to 1).
```{r p1}
P = matrix(c(1/2,1/4,1/4,1/6,1/3,1/2,1/5,1/5,3/5),nrow=3)
colSums(P)
P
```
Here is the corresponding state diagram:
\begin{center}
\begin{tikzpicture}[->,>=stealth,shorten >=1pt,auto,node distance=120pt,thick]
\tikzstyle{every state}=[fill=blue,draw=none,text=white,minimum size=30pt]
\tikzset{every edge/.append style={font=\small}}
\node[state] (A) {$S_1$};
\node[state] (B) [below left of=A] {$S_2$};
\node[state] (C) [below right of=A] {$S_3$};
\path
(A) edge [out=45,in=135,distance=40] node [pos=.5,above] {$0.500$} (A)
(A) edge [bend left=20] node [pos=.4,left] {$0.250$} (B)
(A) edge [bend left=30] node {$0.250$} (C)
(B) edge [out=45+120,in=135+120,distance=40] node [pos=.5,left] {$0.333$} (B)
(B) edge [bend left=30] node {$0.167$} (A)
(B) edge [bend left=20] node [below] {$0.500$} (C)
(C) edge [out=45-120,in=135-120,distance=40] node [pos=.5,right] {$0.600$} (C)
(C) edge [bend left=20] node [pos=.6,right] {$0.200$} (A)
(C) edge [bend left=20] node [pos=.5] {$0.200$} (B);
\end{tikzpicture}
\end{center}
Calculate the eigenstuff:
```{r eigenstuff}
eigenStuff <- eigen(P)
eigenStuff
```
The important information is that, yes, there is a eigenvalue of 1. It is also the largest eigenvalue. The corresponding eigenvector is the first one. We can access easily.
```{r moreeigenstuff}
eigenVecs <- eigenStuff$vectors
steadyStateVec <- eigenVecs[,1]
steadyStateVec
```
This isn't yet a steady state vector for the transition matrix.
We need to **normalize** it so that the sum of the entries is 1.
```{r steadystate}
steadyStateVec <- steadyStateVec/sum(steadyStateVec)
steadyStateVec
```
This is the steady state distribution.
The system will be in $S_1$ with probability `r steadyStateVec[1]`, $S_2$ with probability `r steadyStateVec[2]`, and $S_3$ with probability `r steadyStateVec[3]`.
## Agents again
Let's simulate this Markov chain using _agents_.
We will use a large number of agents and watch how they all move around according to the probabilities defined by the transition matrix $P$.
The idea is simple: if an agent is in state $j$, then it will move probabilistically to its next state according to the probabilities defined by $P[,j]$, the $jth$ column of $P$.
Imagine a large number, $N_1$, of agents in state $S_1$.
We can simulate where they go next using the `sample` function.
Suppose we have $N = 100$ agents in state $S_1$.
```{r next states}
N1 <- 100
nextStates <- sample(1:3, N1, rep=T, prob=P[,1])
nextStates
```
Each of the `r N1` agents ended up in one of $S_1,S_2$ or $S_3$.
Here's how to tally the totals using the `table` function.
Doing so gives us the distribution of next states for these agents.
```{r newDist}
newDist <- as.vector(table(nextStates))
## do this to make sure all three states are represented.
newDist <- c(newDist, rep(0, 3-length(newDist)))
newDist
```
In other words, out of our `r N1` agents,
* `r newDist[1]` ended up in $S_1$,
* `r newDist[2]` ended up in $S_2$,
* `r newDist[3]` ended up in $S_3$,
We could, of course, play the same game with agents in $S_2$ or $S_3$.
Over time, all these agents move between the three states. Letâ€™s simulate this process.
Start with a fixed number of agents arbitrarily distributed between the three states:
```{r numEachState2}
numAgents <- 100000
agents = sample(1:3, numAgents, rep=T)
numEachState <- as.vector(table(agents))
numEachState
```
Hence, initially we have :
* `r numEachState[1]` agents in $S_1$
* `r numEachState[2]` agents in $S_2$
* `r numEachState[3]` agents in $S_3$
Now start them moving around:
```{r nextState123}
nextState1 <- sample(1:3, numEachState[1], rep=T, prob=P[,1])
nextState2 <- sample(1:3, numEachState[2], rep=T, prob=P[,2])
nextState3 <- sample(1:3, numEachState[3], rep=T, prob=P[,3])
```
These aren't much to look at. Looking at the first 50 just shows where these agents ended up.
```{r peek}
head(nextState1, 50)
```
It helps to tally them as before.
```{r fromState123}
fromState1 <- as.vector(table(nextState1))
fromState2 <- as.vector(table(nextState2))
fromState3 <- as.vector(table(nextState3))
```
Each of these represents the *flux* out of each state. For example:
```{r peek fromState1}
fromState1
```
This means that
* `r fromState1[1]` agents went from $S_1$ to $S_1$
* `r fromState1[2]` agents went from $S_1$ to $S_2$
* `r fromState1[3]` agents went from $S_1$ to $S_3$
The other two `fromState` tables contain similar data.
Adding all these together gives the new distribution of the agents after one step of the Markov chain.
```{r numEachState}
numEachState = fromState1 + fromState2 + fromState3
numEachState
```
Hence, after one step of the Markov chain, we have:
* `r numEachState[1]` agents in $S_1$
* `r numEachState[2]` agents in $S_2$
* `r numEachState[3]` agents in $S_3$
Compare this with what we had initially.
Of course, we want to do this a large number of times and see how the agents eventually distribute themselves across the states.
We can do this easily:
```{r compute dist}
## distribute agents uniformly
agents = sample(1:3, numAgents, rep=T)
numEachState <- as.vector(table(agents))
numEachState
numSteps <- 100
for(m in 1:numSteps){
nextState1 <- sample(1:3, numEachState[1], rep=T, prob=P[,1])
nextState2 <- sample(1:3, numEachState[2], rep=T, prob=P[,2])
nextState3 <- sample(1:3, numEachState[3], rep=T, prob=P[,3])
##
fromState1 <- as.vector(table(nextState1))
fromState2 <- as.vector(table(nextState2))
fromState3 <- as.vector(table(nextState3))
##
numEachState <- fromState1 + fromState2 + fromState3
}
numEachState
```
It is better to look at these are proportions and compare to the steady state vector for the transitions matrix $P$.
```{r check SS}
numEachState/numAgents
steadyStateVec
```
Pretty close! In other words, the agents are distributing themselves according the the steady state distribution probabilities. This is as it should be.
Look more closely at all the `fromState` vectors
```{r all fromStates}
fromState1
fromState2
fromState3
```
Notice that the fluxes between the states are nearly balanced---the number of agents leaving an state $S_i$ nearly equals the number coming into $S_i$.
For example, at the end of the simulation, the number of states leaving $S_1$ is given by
```{r flux from 1}
fromState1[2] + fromState1[3]
```
In comparison, the number of agents coming into $S_1$ is
```{r flux to 1}
fromState2[1] + fromState3[1]
```
The numbers would be similar look at $S_2$ or $S_3$. This shows that the system is balanced, indicating that we have reached the steady state distribution of agents.
## Exercises
### 1. Five-State Markov Chain
Consider the following state graph.
\begin{center}
\begin{tikzpicture}[->, >=stealth,shorten < =5pt, shorten > =5pt, auto, thick]
\def \a {72}
\def \n {5}
\def \rad {70pt}
\def \margin {8} % margin in angles, depends on the radius
\tikzstyle{every state}=[fill=blue, draw=none, text=white, minimum size=30pt]
\node[draw, circle,state] (A) at ({360/5*0}:\rad) {$A$};
\node[draw, circle,state] (B) at ({360/5*1}:\rad) {$B$};
\node[draw, circle,state] (C) at ({360/5*2}:\rad) {$C$};
\node[draw, circle,state] (D) at ({360/5*3}:\rad) {$D$};
\node[draw, circle,state] (E) at ({360/5*4}:\rad) {$E$};
\path
(A) edge [bend left] node [pos=0.5] {$0.2$} (B)
(B) edge [bend left] node [pos=0.5] {$0.4$} (C)
(C) edge [bend left] node [pos=0.5] {$0.1$} (D)
(D) edge [bend left] node [pos=0.5] {$0.8$} (E)
(E) edge [bend left] node [pos=0.5] {$0.3$} (A)
(B) edge [bend left] node [pos=0.5] {$0.4$} (A)
(C) edge [bend left] node [pos=0.5] {$0.4$} (B)
(D) edge [bend left] node [pos=0.5] {$0.1$} (C)
(E) edge [bend left] node [pos=0.5] {$0.3$} (D)
(A) edge [bend left] node [pos=0.5] {$0.3$} (E)
(A) edge [out=-45,in=45,distance=60] node [pos=0.5] {$0.5$} (A)
(B) edge [out=-45+\a,in=45+\a,distance=60] node [pos=0.5] {$0.2$} (B)
(C) edge [out=-45+2*\a,in=45+2*\a,distance=60] node [pos=0.5] {$0.5$} (C)
(D) edge [out=-45+3*\a,in=45+3*\a,distance=60] node [pos=0.5] {$0.1$} (D)
(E) edge [out=-45+4*\a,in=45+4*\a,distance=60] node [pos=0.5] {$0.4$} (E);
\end{tikzpicture}
\end{center}
Build the $5 \times 5$ transition matrix for this state graph. Once you have it,
* Find the steady-state distribution.
* Repeat the agents simulation as above. That is, show that, starting with a reasonably large number of agents, that after a few hundred simulations, that the flow of agents in and out of each state are roughly balanced.
### 2. Simulating a Game
Construct a transtion matrix `TransProb` for a game that involves moving through 8 positions, arranged in a circle.
These positions are the states of a Markov chain; we will denote the states as $S_1, S_2, \ldots, S_8$.
$S_1$ in the first position, $S_2$ is the second position, and so on.
In the game, suppose that you roll a dice that tells you how many positions to move.
After you reach $S_8$, the next position is $S_1$, and you go around again.
$$S1 \rightarrow S2 \rightarrow \cdots \rightarrow S8 \rightarrow S1 \rightarrow \cdots$$
Consider two versions of the game:
* **Version 1:** You move from one state to the next rolling a fair *3-sided die*.
Hence if you are in $S_3$ and roll a 2, you move to $S_5$.
If you are in $S_7$ and roll a 3 you move to $S_2$.
* **Version 2:** Same as Version 1, but if you land in $S_5$, then *proceed directly to $S_1$*.
Construct transition matrices for both games.
For each, explore the steady-state behavior.
In other words, what is the probability of being in a particular state, say $S_3$, after a large number of rolls?
Is it different for the the different versions?
Do your work in this notebook, using RMarkdown to explain what you are doing and produce a nice-looking solution in HTML or PDF form.
### 3. Taxi Driver
A city has three different taxi zones, numbered 1, 2, and 3.
A taxi driver operates in all three zones.
The probability that the driver's next passenger has a destination in a particular one of
these zones depends on where the passenger is picked up.
Specifically, whenever the taxi driver is in zone 1, the probability the next passenger is going to zone 1 is .3, to zone 2 is .2, and to zone 3 is .5.
Starting in zone 2, the probability that the next passenger is going to zone 1 is .1, to zone 2 is .8, and to zone 3 is .1.
Finally, whenever the taxi is in zone 3, the probability the next passenger is going to zone 1 is
.4, to zone 2 is .4, and to zone 3 is .2.
Use a Markov chain to simulate the taxi's location in the city, and answer the following questions:
* What are the long-term probabilities of the driver being in each zone of the city?
* If the driver starts the day in zone 3 and transports ten passengers before lunch, what is the probability that the driver will end up in zone 3 for lunch?
* Suppose that the fare for travel between zone 1 and zone 2 is \$5, between zone 1 and zone 3 is \$10, and between zone 2 and zone 3 is \$8. On average, how much money would the taxi driver collect from transporting ten passengers? This may depend on where the driver starts. If so, in which zone should the driver start to maximize the expected revenue from the next ten passengers?