Simulating a Markov Chain in R and Sequence Search

Markov chains are mathematical systems that undergo transitions from one state to another. In this blog post, we will explore how to simulate a Markov chain using R programming language and perform sequence search on the generated data.

Introduction to Markov Chains

A Markov chain is defined as a set of states (S) such that there exists a probability distribution over these states (π), which represents the probability of transitioning from one state to another. The probability of transitioning from state i to state j is represented by P(i,j). We can start at any initial state and then iterate through the states according to their probabilities.

For example, consider a Markov chain with three states: Sunny (S), cloudy (C) and rainy (R). Let’s assume that the transition probabilities are as follows:

  | S | C | R
---------
S | 0.7 | 0.3 | 0.2
C | 0.5 | 0.6 | 0.1
R | 0.2 | 0.2 | 0.8

In this example, the probability of transitioning from Sunny to Cloudy is 0.3 and from Sunny to Rainy is 0.2.

Simulating a Markov Chain in R

To simulate a Markov chain using R programming language, we can use the following code:

## Create a transition matrix for our Markov chain

P = matrix(c(0.7, 0.3, 0.2, 0.2, 0.5, 0.6, 0.1, 0.2, 0.2), 3)

print(P)

This code creates a transition matrix P where each row represents the probabilities of transitioning from one state to another.

Initializing States and Iterating Through Markov Chain

Next, we can initialize our states and iterate through the Markov chain according to their probabilities:

## Initialize the first state

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]

In this code, we first initialize our initial state to “C”. Then we iterate through the Markov chain according to their probabilities and update the current state.

Dividing the Sequence into Groups of Three States

Once we have generated our sequence of states, we can divide it into groups of three states:

## Divide the sequence into groups of three states

collapsed <- paste(states, collapse="")

length(gregexpr("SCC", collapsed)[[1]])

In this code, we use regular expressions to search for “SCC” in our collapsed string. If you want a sliding window (i.e., SCC could occur in position 1-3, or 2-4, etc.), then collapsing the states to a string and to a regex search should do the trick.

Alternatively, 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:

## Divide the sequence into groups of three states

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

sum(chopped == "SCC")

In this code, we use tapply to chop up the sequence into groups of three states and then count the number of those three set states that are equal to SCC.

Conclusion

Simulating a Markov chain using R programming language is a great way to explore complex systems and perform sequence search on generated data. By understanding how to initialize states, iterate through Markov chains according to their probabilities, divide sequences into groups of three states, and use regular expressions or tapply, you can analyze and understand the patterns in your data.

In this blog post, we have covered a basic simulation of a Markov chain with R programming language. However, there are many more ways to simulate complex systems using Markov chains, such as analyzing time series data, modeling network behavior, or understanding population dynamics.

Last modified on 2023-07-10