---
title: "Introduction to Markov Chains"
output: html_document
---
# Review of Matrix Manipulations
## Getting Started
Define a couple matrices
```{r matrices1}
mat1 <- matrix(c(3,2,-1,4), nrow=2)
mat2 <- matrix(c(0,1,2,-1), nrow=2)
mat1
mat2
```
## Multiplication
Look at them and multiply together. Note the %*% means matrix multiplication
```{r matrices 2}
mat1
mat2
mat1 %*% mat2
```
## Vectors
Vectors work the same way, the are just matrices with one column.
```{r vectors 1}
vec1 <- matrix(c(3,-2), nrow=2)
vec1
```
Take a peek at what we have.
```{r vectors 2}
mat1
vec1
mat1 %*% vec1
```
## Determinants, inverses
You can take determinants, inverses, etc.
```{r determinants}
det(mat1)
mat1_inv <- solve(mat1)
mat1 %*% mat1_inv
```
**Question:**
*Why does `solve` return the inverse of matrix A?*
What is being *solved* here?
Try the following:
```{r}
sol <- solve(mat1, vec1)
```
What did this function do?
Does this shed any light on why `solve` would return the inverse of a matrix?
# Markov Chains
Here we do a quick tour of Markov chains (assuming you've already seen Markov chains in Linear Algebra).
## Weather in the Land of Oz.
There are only three possible states:
* Rain
* Nice
* Snow
The weather changes every day, with the following probabilities:
* If it's Nice, then it will either Rain or Snow the next day with equal probability.
* If it's Rain, then the next day it will Rain again with probability 1/2, otherwise it will be either Nice or Snow with equal probability.
* If it's Snow, then the next day it will Snow again with probability 1/2, otherwise it will be either Nice or Rain with equal probability.
From this we define the **Transition Matrix** $P$.
The entry in row $i$, column $j$ gives the probability of transitioning from state $j$ to state $i$.
We will use labels R (Rain), N (Nice), and S (Snow) for the rows and columns.
This isn't necessary but it does make it easier to see what is going on.
```{r transition 1 }
TransProb <- matrix(c(1/2,1/4,1/4,1/2,0,1/2,1/4,1/4,1/2),nrow=3)
rownames(TransProb) <- c("R","N","S")
colnames(TransProb) <- c("R","N","S")
TransProb
```
*Note that each column sums to 1, but not each row.*
**Question**: In the long run, what is the probability of each type of weather for a given day?
## Approach 1: Steady-state distribution
Suppose on Day 1, it is R raining, what is the probability of the weather the next day?
Answer:
* Prob(Rain) = 0.5
* Prob(Nice) = 0.25
* Prob(Snow) = 0.25
One way to get this is via matrix multiplication:
```{r day vectors}
Day0Vec <- matrix(c(1,0,0), nrow=3)
Day0Vec
Day1Vec <- TransProb %*% Day0Vec
Day1Vec
```
We can continue this for several days.
*Note: the parentheses around the assignment means print the assignment values.*
```{r day vectors 2}
(Day2Vec <- TransProb %*% Day1Vec)
(Day3Vec <- TransProb %*% Day2Vec)
```
Or a hundred days.
```{r day vectors 100 days}
NumDays <- 100
DayVec <- Day0Vec
print(DayVec)
for(i in 1:NumDays){
DayVec <- TransProb %*% DayVec
}
DayVec
```
Or, use a different starting day.
Let's start with a random probability distribution on the three states.
```{r day vector new start}
DayVec <- c(1,3,5)
#make sure the probabilities add up to one
DayVec <- DayVec/sum(DayVec)
print(DayVec)
for(i in 1:NumDays){
DayVec <- TransProb %*% DayVec
}
DayVec
```
We have found the *steady-state* distribution.
In the long run, the weather stabilizes to 40% R, 20% N, and 40% S.
You can check, not matter what the starting distribution is, you end up in the same steady-state.
## Approach 2: Matrix Powers
Since this involves powers of TransProb matrix, let's look at these.
R doesn't have a matrix exponentiation command (there is a package that includes one, namely **expm**).
We can compute powers the direct way, by repeated multiplication.
Here is `TransProb^25`.
```{r transition powers}
TransProb2 <- TransProb
M <- 25
print(TransProb2)
for(i in 1:M){
TransProb2 <- TransProb2 %*% TransProb
}
TransProb2
```
There those numbers are again!
## Approach 3: Agents
Another way to look at this, follow the weather from day to day.
Day 1: Snow = S.
Day 2: R with probability 1/4, N with probability 1/4, S with probability 1/2
Let's determine the next day's weather using these probabilities.
Use the `sample` function, with weights given by the probabilities.
```{r agents 1}
weather <- c("R","N","S")
sample(weather, 1, prob=c(1/4,1/4,1/2))
w0 <- "S"
probs <- TransProb[,w0]
sample(weather, 1, prob=probs)
```
Now do this over and over.
Start again in Snow, and use the probabilities to determine the next day's weather.
```{r agents 2}
##the initial weather state
currWeather <- "S"
probs <- TransProb[,currWeather]
probs
##the next weather state
currWeather <- sample(weather,1,prob=probs)
currWeather
```
Now use the current weather, along with the probabilites, to determine the next day's weather.
```{r agents 3}
currWeather
probs <- TransProb[,currWeather]
probs
##the next weather state
currWeather <- sample(weather,1,prob=probs)
currWeather
```
Repeat....
```{r agents 4}
currWeather
probs<-TransProb[,currWeather]
probs
##the next weather state
currWeather <- sample(weather,1,prob=probs)
currWeather
```
Repeat....wait, let's use a loop.
Start with Snow, then simulate 10000 days.
Use each day's weather, together with the probabilites, to determine the next day's weather.
```{r agents loop}
numDays <- 10000
currWeather <- "S"
all_weather <- c()
for(i in 1:numDays){
all_weather[i] <- currWeather
probs <- TransProb[,currWeather]
currWeather <- sample(weather,1,prob=probs)
}
table(all_weather)/numDays
```
There are those numbers again!
## Approach 4: Eigenvectors and Eigenvalues
If $v$ is a steady state vector of matrix `TransProb`, then `TransProb`$\cdot v=v$.
In other words, $v$ is an *eigenvector* of matrix `TransProb` with *eigenvalue* 1.
**Fact**: A well-behaved transition matrix always has a steady-state vector, hence it has a eigenvector with eigenvalue 1.
Here's how to get eigen state information in R:
```{r eigenvectors}
eState <- eigen(TransProb)
eState
```
We can pull out the eigenvalues and vectors
```{r eigenvectors values}
(eVals <- eState$value)
(eVecs <- eState$vectors)
```
Note the first eigenvalue is 1, as predicted. The corresponding eigenvalue is
```{r eigenvector values 2}
(eVec1 <- eVecs[,1])
```
Eigenvectors are only unique up to scale. In this case, the vector needs to be rescaled to be a probability distribution---i.e., all non-negative entries that sum to 1.
```{r eigenvector values 3}
(eVec1 <- eVec1/sum(eVec1))
```
There are those numbers again!
# Exercise 1
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 S1, S2, ..., S8.
S1 in the first position, S2 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 S8, the next position is S1, and you go around again.
S1->S2-> ... ->S8->S1->etc
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 S3 and roll a 2, you move to S5.
If you are in S7 and roll a 3 you move to S2.
* **Version 2:** Same as Version 1, but if you land in S5, then *proceed directly to S1*.
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 S3, 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.
# Exercise 2
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:
1. What are the long-term probabilities of the driver being in each zone of the city?
2. 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?
3. 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?