--- 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?