Generating Markov transition matrix in Python

Question

Imagine I have a series of 4 possible Markovian states (A, B, C, D):

X = [A, B, B, C, B, A, D, D, A, B, A, D, ....]

How can I generate a Markov transformation matrix using Python? The matrix must be 4 by 4, showing the probability of moving from each state to the other 3 states. I've been looking at many examples online but in all of them, the matrix is given, not calculated based on data. I also looked into hmmlearn but nowhere I read on how to have it spit out the transition matrix. Is there a library that I can use for this purpose?

Here is an R code for the exact thing I am trying to do in Python: https://stats.stackexchange.com/questions/26722/calculate-transition-matrix-markov-in-r

In raw Python, you would need to use a list of lists. This sort of thing is more naturally done in `numpy` or `pandas`. If you want to use one of those tools, perhaps you can add the appropriate tag. In any event, what is the input of your problem? A finite list of states? — John Coleman, Oct 10 '17 at 01:18
"the probability of moving from each state to the other 3 states" Shouldn't it rather be "to any of the 4 states", since a state can persist for one or more time steps? — NerdOnTour, Jan 03 '22 at 08:56

John Coleman · Accepted Answer · 2017-10-12T15:46:58.403

This might give you some ideas:

transitions = ['A', 'B', 'B', 'C', 'B', 'A', 'D', 'D', 'A', 'B', 'A', 'D']

def rank(c):
    return ord(c) - ord('A')

T = [rank(c) for c in transitions]

#create matrix of zeros

M = [[0]*4 for _ in range(4)]

for (i,j) in zip(T,T[1:]):
    M[i][j] += 1

#now convert to probabilities:
for row in M:
    n = sum(row)
    if n > 0:
        row[:] = [f/sum(row) for f in row]

#print M:

for row in M:
    print(row)

output:

[0.0, 0.5, 0.0, 0.5]
[0.5, 0.25, 0.25, 0.0]
[0.0, 1.0, 0.0, 0.0]
[0.5, 0.0, 0.0, 0.5]

On Edit Here is a function which implements the above ideas:

#the following code takes a list such as
#[1,1,2,6,8,5,5,7,8,8,1,1,4,5,5,0,0,0,1,1,4,4,5,1,3,3,4,5,4,1,1]
#with states labeled as successive integers starting with 0
#and returns a transition matrix, M,
#where M[i][j] is the probability of transitioning from i to j

def transition_matrix(transitions):
    n = 1+ max(transitions) #number of states

    M = [[0]*n for _ in range(n)]

    for (i,j) in zip(transitions,transitions[1:]):
        M[i][j] += 1

    #now convert to probabilities:
    for row in M:
        s = sum(row)
        if s > 0:
            row[:] = [f/s for f in row]
    return M

#test:

t = [1,1,2,6,8,5,5,7,8,8,1,1,4,5,5,0,0,0,1,1,4,4,5,1,3,3,4,5,4,1,1]
m = transition_matrix(t)
for row in m: print(' '.join('{0:.2f}'.format(x) for x in row))

Output:

0.67 0.33 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.50 0.12 0.12 0.25 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00
0.00 0.00 0.00 0.50 0.50 0.00 0.00 0.00 0.00
0.00 0.20 0.00 0.00 0.20 0.60 0.00 0.00 0.00
0.17 0.17 0.00 0.00 0.17 0.33 0.00 0.17 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00
0.00 0.33 0.00 0.00 0.00 0.33 0.00 0.00 0.33

Thanks, John. The only thing is that the transition matrix resulting from your code is not exactly the Markov transition matrix. Your code is dividing the row cells by n=11. Is there anyway you can conveniently fix the code to give the table in the screenshot I added to the question's **"update"** section? — st19297, Oct 10 '17 at 16:23
@st19297 I just replaced the global `n` with row-specific `n` (making the entries conditional probabilities). Since I don't like to divide by 0, the above code leaves a row of zeros unchanged. It is impossible to estimate transition probabilities from a given state when no transitions from that state have been observed. — John Coleman, Oct 10 '17 at 17:08
@JohnColeman, I have updated your code using numpy (it's shorter) https://stackoverflow.com/a/72673234/7564809 but I haven't checked for zero division. As long as I know each row in a transition matrix should add up to one, shouldn't it? — Andrea Dalseno, Jun 18 '22 at 22:49

Iain D · Answer 2 · 2020-09-29T15:30:06.767

If you want to do it all in pandas, here is an approach that works for non numeric data:

import pandas as pd
transitions = ['A', 'B', 'B', 'C', 'B', 'A', 'D', 'D', 'A', 'B', 'A', 'D']

df = pd.DataFrame(transitions)

# create a new column with data shifted one space
df['shift'] = df[0].shift(-1)

# add a count column (for group by function)
df['count'] = 1

# groupby and then unstack, fill the zeros
trans_mat = df.groupby([0, 'shift']).count().unstack().fillna(0)

# normalise by occurences and save values to get transition matrix
trans_mat = trans_mat.div(trans_mat.sum(axis=1), axis=0).values

It's slower than the pure python approach but maybe worth it for flexibility and to avoid creating your own function.

Great solution, thanks. For those interested in the order of the transitions in their matrix: follow the order of `trans_mat` before it is being normalized. — sander, Jun 29 '23 at 13:49

Polykarpos Polykarpidis · Answer 3 · 2021-07-09T10:37:39.047

The following code provides another solution about Markov transition matrix order 1. Your data can be list of integers, list of strings, or a string. The negative think is that this solution -most likely- requires time and memory.

creates a Markov transition matrix order 1 (bigrams)
generates 1000 integers in order to train the Markov transition matrix to a dataset.
train the Markov transition matrix

Until here we have the solution of the question. The following code try to solve an additional problem. Specifically, the generating data according to the trained Markov task.

transform probabilities of markov transition matrix to cumulative (arithmetic coding)
generating 30 data

import pandas as pd

def transition_matrix_order1(data):
    alphabet = []
    for element in data:
        if element not in alphabet:
            alphabet.append(element)
    alphabet.sort()
    
    previous = data[0]
    matrix = pd.DataFrame(0.0, index=alphabet, columns=alphabet)
    
    for i in data[1:]:
        matrix[i][previous]    += 1.0
        previous = i
    
    total = matrix.sum()
    for element in alphabet:
        matrix[element] = matrix.div(total[element])[element]
    
    return matrix, alphabet



#create data using random integers========
import random
data = [random.randint(1,5) for i in range(1000)] #You can also put list of strings or a string as input data



#create markov transition matrix order 1 (bigram)
markov_matrix, alphabet = transition_matrix_order1(data)



#=the following code uses the probabilities in order to create new data.=



#transform probabilities of markov transition matrix to cumulative
for column in alphabet:
    for pos, index in enumerate(alphabet[1:]):
        markov_matrix[column][index] += markov_matrix[column][alphabet[pos]]




#generating 30 data
generated_data = []
feed = random.choice(alphabet)
generated_data.append(feed)
for i in range(30):
    random_value = random.uniform(0, 1)
    for i in alphabet:
        if markov_matrix[feed][i] >= random_value:
            generated_data.append(i)
            feed = i
            break



print(generated_data)

To anyone willing to use this solution please beware that Pandas' DataFrame is inverted regarding a normal matrix. In other words, the variable `matrix` returned by the `transition_matrix_order1()` function in the code above is not a functional transition matrix. In order to obtain a real transition matrix, transposing the `matrix` variable is enough. That is, `matrix = = matrix.transpose()`. — Razvi, Oct 25 '22 at 21:34

sander · Answer 4 · 2023-07-20T09:08:03.870

In Pandas there is a much easier solution: pd.crosstab. Given your sequence:

X = ["A", "B", "B", "C", "B", "A", "D", "D", "A", "B", "A", "D"]

matrix = pd.crosstab(
    pd.Series(X[:-1], name='from'),
    pd.Series(X[1:], name='to'),
    normalize=0
)

Resulting in the following pd.DataFrame:

    to  A   B    C    D
from                
A       0.0 0.50 0.00 0.5
B       0.5 0.25 0.25 0.0
C       0.0 1.00 0.00 0.0
D       0.5 0.00 0.00 0.5

If you want a np.array instead, use matrix.to_numpy() which results in:

[[0.   0.5  0.   0.5 ]
 [0.5  0.25 0.25 0.  ]
 [0.   1.   0.   0.  ]
 [0.5  0.   0.   0.5 ]]

score 0 · Answer 5 · answered Jun 18 '22 at 22:43

Thank you @john-coleman , I have updated your code using numpy:

import numpy as np

def transition_matrix(transitions):
    n = 1+ max(transitions) #number of states

    M = np.zeros((n,n))

    for (i,j) in zip(transitions,transitions[1:]):
        M[i][j] += 1

    #now convert to probabilities:
    M = M/M.sum(axis=1, keepdims=True)
    return M

t = [1,1,2,6,8,5,5,7,8,8,1,1,4,5,5,0,0,0,1,1,4,4,5,1,3,3,4,5,4,1,1]
m = transition_matrix(t)
for row in m: print(' '.join(f'{x:.2f}' for x in row))

The output is the same:

0.67 0.33 0.00 0.00 0.00 0.00 0.00 0.00 0.00
0.00 0.50 0.12 0.12 0.25 0.00 0.00 0.00 0.00
0.00 0.00 0.00 0.00 0.00 0.00 1.00 0.00 0.00
0.00 0.00 0.00 0.50 0.50 0.00 0.00 0.00 0.00
0.00 0.20 0.00 0.00 0.20 0.60 0.00 0.00 0.00
0.17 0.17 0.00 0.00 0.17 0.33 0.00 0.17 0.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00
0.00 0.00 0.00 0.00 0.00 0.00 0.00 0.00 1.00
0.00 0.33 0.00 0.00 0.00 0.33 0.00 0.00 0.33

Generating Markov transition matrix in Python

5 Answers5

Linked