Consider a markov chain of 4 states $\{S_1, S_2, S_3, S_4\}$ described by the transition matrix

$$ A = \begin{bmatrix} .25 & .20 & .25 & .30 \\ .20 & .30 & .25 & .30 \\ .25 & .20 & .40 & .10 \\ .30 & .30 & .10 & .30 \end{bmatrix} $$

I understand that given an initial state $x^{(0)}$ the probability that the system is in state $i$ after $n$ steps is

$$ x^{(n)} = A^n x^{(0)} $$

However, I want to know how to generate a realization of a sequence states. For example, from the matrix above one possible sequence of states over $n = 0, 1, 2, 3$ steps could be

$$ S_1 \rightarrow S_3 \rightarrow S_4 \rightarrow S_1 $$

Is there a way to programmatically generate a particular sequence or realization from the transition matrix, either based on a random initial state or given an initial state $x^{(0)}$?

  • $\begingroup$ In matlab you can use the functions dtmc and simulate. $\endgroup$ – L Y Mar 30 at 2:33
  • $\begingroup$ Thanks but I'm looking for a general algorithm or approach, not just a function call. $\endgroup$ – mikeck Mar 30 at 2:38
  • $\begingroup$ Although implementation is often mixed with substantive content in questions, we are supposed to be a site for providing information about statistics, machine learning, etc., not code. It can be good to provide code as well, but please elaborate your substantive answer in text for people who don't read this language well enough to recognize & extract the answer from the code. $\endgroup$ – gung - Reinstate Monica Mar 30 at 4:32

Starting from state $i$, draw a sample according to the distribution defined by that column and update the state according to the output.

Here's an implementation in Python using $0$-index.

import numpy as np

A = np.array([[.25, .2, .25, .3], [.2, .3, .25, .3], [.25, .2, .4, .1], [.3, .3, .1, .3]])
xk = np.arange(len(A))

def generate_sample(cur_state):
    return np.random.choice(xk, 1,  p = A[:, cur_state])

initial_state = 0
sample_len = 100
output = [-1 for i in range(sample_len)]
output[0] = initial_state
for i in range(1, sample_len):
    output[i] = generate_sample(output[i-1])[0]

Remark: In most literature, each row sums to $1$ but for your example, each column sums to $1$.

| cite | improve this answer | |
  • $\begingroup$ Thanks, I started thinking of an approach like this after I posted the question. I was hoping there would be a way to do this without a for loop but it's probably necessary to do it like this since we need to track the memory of states. I'll be doing this in R but will be using massive (sparse) matrices, so I'll probably use C++ via Rcpp for the implementation. $\endgroup$ – mikeck Mar 30 at 15:39

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.