generate a realization from a transition matrix

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)}$$?

• In matlab you can use the functions dtmc and simulate.
– L Y
Commented Mar 30, 2020 at 2:33
• Thanks but I'm looking for a general algorithm or approach, not just a function call. Commented Mar 30, 2020 at 2:38
• 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. Commented Mar 30, 2020 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]
print(output)


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

• 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. Commented Mar 30, 2020 at 15:39