When working with Markov chains, it is often necessary to calculate the stationary distribution. The stationary distribution represents the long-term probabilities of being in each state of the Markov chain. In this article, we will explore three different ways to calculate the stationary distribution of a Markov chain in Python.

Option 1: Using NumPy

One way to calculate the stationary distribution is by using the NumPy library. NumPy provides various mathematical functions and tools that are useful for working with arrays and matrices. Here is a sample code that demonstrates how to calculate the stationary distribution using NumPy:

import numpy as np

# Define the transition matrix
transition_matrix = np.array([[0.2, 0.8],
                              [0.6, 0.4]])

# Calculate the eigenvalues and eigenvectors of the transpose of the transition matrix
eigenvalues, eigenvectors = np.linalg.eig(transition_matrix.T)

# Find the eigenvector corresponding to the eigenvalue of 1
stationary_distribution = np.real(eigenvectors[:, np.isclose(eigenvalues, 1)])

# Normalize the stationary distribution
stationary_distribution /= np.sum(stationary_distribution)

print(stationary_distribution)

In this code, we first define the transition matrix of the Markov chain. Then, we use the NumPy function np.linalg.eig to calculate the eigenvalues and eigenvectors of the transpose of the transition matrix. We find the eigenvector corresponding to the eigenvalue of 1, which represents the stationary distribution. Finally, we normalize the stationary distribution by dividing it by the sum of its elements.

Option 2: Using SciPy

Another way to calculate the stationary distribution is by using the SciPy library. SciPy is a powerful library for scientific computing in Python, and it provides various functions for linear algebra, optimization, and more. Here is a sample code that demonstrates how to calculate the stationary distribution using SciPy:

import numpy as np
from scipy.linalg import eig

# Define the transition matrix
transition_matrix = np.array([[0.2, 0.8],
                              [0.6, 0.4]])

# Calculate the eigenvalues and eigenvectors of the transpose of the transition matrix
eigenvalues, eigenvectors = eig(transition_matrix.T)

# Find the eigenvector corresponding to the eigenvalue of 1
stationary_distribution = np.real(eigenvectors[:, np.isclose(eigenvalues, 1)])

# Normalize the stationary distribution
stationary_distribution /= np.sum(stationary_distribution)

print(stationary_distribution)

In this code, we import the eig function from the scipy.linalg module to calculate the eigenvalues and eigenvectors of the transpose of the transition matrix. The rest of the code is similar to Option 1.

Option 3: Using SymPy

A third way to calculate the stationary distribution is by using the SymPy library. SymPy is a Python library for symbolic mathematics, and it provides various functions for algebraic manipulations, solving equations, and more. Here is a sample code that demonstrates how to calculate the stationary distribution using SymPy:

import sympy as sp

# Define the transition matrix
transition_matrix = sp.Matrix([[0.2, 0.8],
                              [0.6, 0.4]])

# Calculate the eigenvalues and eigenvectors of the transpose of the transition matrix
eigenvalues, eigenvectors = transition_matrix.T.eigenvects()

# Find the eigenvector corresponding to the eigenvalue of 1
stationary_distribution = np.real(eigenvectors[0][2][0])

# Normalize the stationary distribution
stationary_distribution /= sp.Sum(stationary_distribution).doit()

print(stationary_distribution)

In this code, we first import the sympy module as sp. We define the transition matrix as a SymPy Matrix object. Then, we use the eigenvects method to calculate the eigenvalues and eigenvectors of the transpose of the transition matrix. The rest of the code is similar to Option 1 and 2.

After exploring these three options, it is clear that Option 1, which uses NumPy, is the most straightforward and efficient way to calculate the stationary distribution of a Markov chain in Python. NumPy provides a wide range of mathematical functions and tools that are specifically designed for numerical computations, making it the ideal choice for this task.