In [3]:
import numpy as np
import random

$k$-step transition probability¶

1. by sampling¶

In [4]:
def update_state(s):
    """
        Does one step of the Markov chain. 
    """
    
    if s==0:
        if random.random() < 0.99: return 0
        else: return 1
    elif s==1:
        if random.random() < 0.9: return 0
        else: return 1
In [5]:
def many_updates(s,k):
    for _ in range(k):
        s = update_state(s)
        
    return s

First, we sample $k$ steps of the Markov chain. We do this $10^5$ times and thus sample the probability distribution at the $k$-th step. We compare two case: one where we start with $s=0$ and another starting with $s=1$.

In [6]:
n = 10**5
L = [many_updates(s=0,k=100) for _ in range(n)]

print(L.count(0)/n, L.count(1)/n)

L = [many_updates(s=1,k=100) for _ in range(n)]

print(L.count(0)/n, L.count(1)/n)

0.98883 0.01117
0.98904 0.01096

2. By Chapman-Kolmogorov formula¶

Next, we compute the probability exactly using formula from class: by computing powers of the matrix.

In [7]:
P = np.matrix([[0.99, .01], [.9, .1]]); P
Out[7]:
matrix([[0.99, 0.01],
        [0.9 , 0.1 ]])
In [8]:
P*P
Out[8]:
matrix([[0.9891, 0.0109],
        [0.981 , 0.019 ]])
In [9]:
P**100
Out[9]:
matrix([[0.98901099, 0.01098901],
        [0.98901099, 0.01098901]])
In [ ]: