CODEWITHSUNDEEP
rprobabilitymarkov-chains
Kyle Knudson
Kyle Knudson

Reputation: 5

Simulating a Markov Chain in R and Sequence Search

So I am working on simulating a Markov Chain with R in which the states are Sunny (S), cloudy (C) and rainy (R) and am looking to figure out the probability that a sunny day is followed by two consecutive cloudy days.

Here is what I have so far:

    P = matrix(c(0.7, 0.3, 0.2, 0.2, 0.5, 0.6, 0.1, 0.2, 0.2), 3)
    print(P)
    x = c("S", "C", "R") 
    n = 10000
    states = character(n+100)
    states[1] = "C"

    for (i in 2:(n+100)){
    if (states[i-1] == "S") {cond.prob = P[1,]}
    else if (states[i-1] == "C") {cond.prob = P[2,]}
    else {cond.prob = P[3,]}
    states[i]=sample(x, 1, prob = cond.prob )
    }

    print(states[1:100])
    states = states[-(1:100)]  
    head(states)
    tail(states)
    states[1:200]

At the end of this I am left with a sequence of states. I am looking to divide this sequence into groups of three states (For the three days in the chain) and then count the number of those three set states that are equal to SCC.

I am running a blank on how I would go about doing this and any help would be greatly appreciated!!

Upvotes: 0

Views: 233

Answers (2)

David Heckmann
David Heckmann

Reputation: 2939

Eric provided the correct answer, but just for the sake of completeness: you can get the probability you are looking for using the equilibrium distribution of the chain:

# the equilibrium distribution
e <- eigen(t(P))$vectors[,1]
e.norm <- e/sum(e)
# probability of SCC
e.norm[1] * P[1,2] * P[2,2]

This is less computationally expensive and will give you a more accurate estimate of the probability, because your simulation will be biased towards the initial state of the chain.

Upvotes: 1

Eric
Eric

Reputation: 1701

Assuming you want a sliding window (i.e., SCC could occur in position 1-3, or 2-4, etc.), collapsing the states to a string and to a regex search should do the trick:

collapsed <- paste(states, collapse="")
length(gregexpr("SCC", collapsed)[[1]])

On the other hand, if you DON'T want a sliding window (i.e., SCC must be in position 1-3, or 4-6, or 7-9, etc.) then you can chop up the sequence using tapply:

indexer <- rep(1:(ceiling(length(states)/3)), each=3, length=length(states))
chopped <- tapply(states, indexer, paste0, collapse="")
sum(chopped == "SCC")

Upvotes: 1

Related Questions