CODEWITHSUNDEEP
rmarkov-chains
G. Debailly
G. Debailly

Reputation: 1331

Markov chain in R not leaving initial state

My Markov chain simulation will not leave the initial state 1. The 4x4 transition matrix has absorption states 0 and 3.

The same code is working for a 3x3 transition matrix without absorption states. What is wrong?

A plot from my code: Plot of Markov chain simulation. And hopefully a working example:

# Building Markov transition matrix)
Pc <- matrix(c(1, 0.05, 0, 0, 0, 0.85, 0.05, 0, 0, 0.1, 0.65, 0, 0, 0, 0.3, 1), 
            nrow = 4, ncol = 4)
rownames(Pc) <- c(0,1,2,3)
colnames(Pc) <- c(0,1,2,3)

#Simulating Markov chain 
markovSimulation <- function(matrix, length, initialState) {
    chain <- array(0,c(length,1))
    chain[1] <- initialState
    for (j in 2:length){
        chain[j] <- sample(1:4, size = 1, prob=matrix[chain[length(chain)] + 1, ])
    }#for loop
    return(chain)
}#markovSimulation

# Calling simulator and plotting result
simulatedChain <- markovSimulation(Pc, 10, 1)
plot(simulatedChain)

Upvotes: 2

Views: 660

Answers (2)

IRTFM
IRTFM

Reputation: 263342

I saw the confusion as due to the fact that you created character indices for your transition matrix and then didn't use them as intended. Using 1 as an i-index would give the first row. Using "1" would give the second row. This uses a character vector of states and uses the character value as a row index into your Pc matrix. Better is to use tmat inside the function rather than a name matrix since that is also a function name:

markovSimulation <- function(tmat, length, initialState) {
     chain <- array(0,c(length,1))
     chain[1] <- initialState
     for (j in 2:length){
         chain[j] <- sample(0:3, size = 1, prob=tmat[ chain[j-1], ] )
     }#for loop
     return(chain)
 }#markovSimulation

 # Calling simulator and plotting result
 simulatedChain <- markovSimulation(Pc, 100, "1")
 plot(simulatedChain)

Upvotes: 1

Jan Sila
Jan Sila

Reputation: 1593

welcome to SO! Good effort, but I see you got two things wrong:

chain[j] <- sample(1:4, size = 1, prob=matrix[chain[length(chain)] + 1, ])

you need to sample from c(0,1,2,3) as you index your states from 0. and the probabilities need to depend on the previous state, so prob=matrix[chain[j-1]+1].

Then I'm getting reasonable results like this: 

    > markovSimulation(Pc,10,1)
      [,1]
 [1,]    1
 [2,]    1
 [3,]    2
 [4,]    3
 [5,]    3
 [6,]    3
 [7,]    3
 [8,]    3
 [9,]    3
[10,]    3

+1 for nice first post, good reproducible example!

Upvotes: 1

Related Questions