In [1]:
# Code attribution: Yiyin Shen, Tyler Caraza-Harter
# Imports
import numpy
In [2]:
# Simulate the Markov chain that models the weather [sun, cloud, rain]
rng = numpy.random.default_rng(seed = 42)
tr = numpy.array([[0.5, 0.3, 0.2], [0.2, 0.5, 0.3], [0.3, 0.2, 0.5]])
seq = [rng.choice(3)]
for i in range(1000):
seq.append(rng.choice(3, p = tr[seq[i], :]))
seq[:20]
Out[2]:
[0, 0, 2, 2, 0, 2, 2, 2, 0, 0, 0, 2, 2, 2, 1, 1, 1, 0, 2, 2]
In [3]:
# Estimate the transitions again
tr_hat = numpy.array([[0] * 3] * 3)
for i in range(999):
tr_hat[seq[i], seq[i + 1]] += 1
tr_hat = tr_hat / numpy.sum(tr_hat, axis = 1).reshape(-1, 1)
tr_hat
Out[3]:
array([[0.51714286, 0.25428571, 0.22857143],
[0.22077922, 0.49350649, 0.28571429],
[0.29618768, 0.19648094, 0.50733138]])
In [4]:
# Estimate the stationary distribution
mu_hat, _ = numpy.histogram(seq, bins = 3)
mu_hat / numpy.sum(mu_hat)
Out[4]:
array([0.35164835, 0.30769231, 0.34065934])
In [5]:
# Approximate the stationary distribution: anything * transition * transition * transition ...
mu_hat = [1, 0, 0]
for i in range(100):
mu_hat = mu_hat @ tr
mu_hat
Out[5]:
array([0.33333333, 0.33333333, 0.33333333])
In [6]:
# Verify the stationary distribution: mu * transition = mu
[1/3] * 3 @ tr
Out[6]:
array([0.33333333, 0.33333333, 0.33333333])