Markov Chains



In [1]:

    
import numpy as np



In [2]:

    
def stationary_distribution(P: np.array) -> np.array:
    A = P - np.eye(len(P))
    A[:, (len(P) - 1)] = np.ones([len(P)])
    
    p_stationary = np.linalg.inv(A)[len(P) - 1, :]
    return p_stationary



In [3]:

    
def check_detailed_balance(transition, pi=None):
    if pi is None:
        pi = stationary_distribution(transition)
    detailed_balance= True
    for x in range(4):
        for y in range(4):
            detailed_balance &= (pi[x] * transition[x, y] == pi[y] * transition[y, x])
    return detailed_balance

Exercises Durret Book

Exercise 1.2



In [4]:

    
transition_2 = np.array([[     0,     1,     0,     0,     0,     0],
                         [  1/25,  8/25, 16/25,     0,     0,     0],
                         [     0,  4/25, 12/25,  9/25,     0,     0],
                         [     0,     0,  9/25, 12/25,  4/25,     0],
                         [     0,     0,     0, 16/25,  8/25,  1/25],
                         [     0,     0,     0,     0,     1,     0]])

Exercise 1.3



In [5]:

    
transition_3 = (1/16) * np.array([[ 3, 2, 2, 2, 3, 4],
                                  [ 4, 3, 2, 2, 2, 3],
                                  [ 3, 4, 3, 2, 2, 2],
                                  [ 2, 3, 4, 3, 2, 2],
                                  [ 2, 2, 3, 4, 3, 2],
                                  [ 2, 2, 2, 3, 4, 3]])

Exercise 1.7



In [6]:

    
transition_7 = np.array([[0.6, 0.4,   0,   0],
                         [  0,   0, 0.6, 0.4],
                         [0.6, 0.4,   0,   0],
                         [  0,   0, 0.3, 0.7]])

Exercise 1.9



In [ ]:

Exercise 1.10



In [7]:

    
transition_10c = np.array([[0.7, 0, 0.3, 0], 
                           [0.2, 0.5, 0.3, 0], 
                           [0.1, 0.2, 0.1, 0.3], 
                           [0, 0.4, 0, 0.6]])
print(stationary_distribution(transition_10c))









    



[0.25454545 0.29090909 0.18181818 0.27272727]

Exercise 1.11



In [8]:

    
print(stationary_distribution(transition_2))









    



[0.00396825 0.09920635 0.3968254  0.3968254  0.09920635 0.00396825]



In [9]:

    
print(stationary_distribution(transition_3))









    



[0.16666667 0.16666667 0.16666667 0.16666667 0.16666667 0.16666667]



In [10]:

    
print(stationary_distribution(transition_7))









    



[0.31034483 0.20689655 0.20689655 0.27586207]

Exercise 1.12



In [11]:

    
transition_12 = np.array([[  0, 2/3,   0, 1/3],
                          [1/3,   0, 2/3,   0],
                          [  0, 1/6,   0, 5/6],
                          [2/5,   0, 3/5,   0]])



In [12]:

    
print(stationary_distribution(transition_12))









    



[0.18817204 0.17741935 0.31182796 0.32258065]

Exercise 1.13



In [13]:

    
transition_13 = np.array([[0.0, 0.0, 0.1, 0.9],
                          [0.0, 0.0, 0.6, 0.4],
                          [0.8, 0.2, 0.0, 0.0],
                          [0.4, 0.6, 0.0, 0.0]])



In [14]:

    
stationary_distribution(np.linalg.matrix_power(transition_13, 3))









    Out[14]:





array([0.26666667, 0.23333333, 0.16666667, 0.33333333])

Exercise: Ruiz Family



In [15]:

    
transition_ruiz = np.array([[0.0, 1.0, 0.0, 0.0, 0.0],
                            [0.3, 0.0, 0.7, 0.0, 0.0],
                            [0.3, 0.0, 0.0, 0.7, 0.0],
                            [0.3, 0.0, 0.0, 0.0, 0.7],
                            [1.0, 0.0, 0.0, 0.0, 0.0]])



In [16]:

    
np.linalg.matrix_power(transition_ruiz, 4)









    Out[16]:





array([[0.237, 0.21 , 0.21 , 0.   , 0.343],
       [0.469, 0.237, 0.147, 0.147, 0.   ],
       [0.126, 0.58 , 0.147, 0.147, 0.   ],
       [0.273, 0.09 , 0.49 , 0.147, 0.   ],
       [0.21 , 0.3  , 0.   , 0.49 , 0.   ]])



In [17]:

    
stationary_distribution(transition_ruiz)









    Out[17]:





array([0.28304557, 0.28304557, 0.1981319 , 0.13869233, 0.09708463])



In [ ]:

Exercise 1.19



In [18]:

    
transition_19 = np.array([[0.7, 0.3],
                          [0.1, 0.9]])

Exercise 1.19 a)



In [19]:

    
transition_19_3 = np.linalg.matrix_power(transition_19, 3)
transition_19_3









    Out[19]:





array([[0.412, 0.588],
       [0.196, 0.804]])

Exercise 1.19 b)



In [20]:

    
np.array([[0.75, 0.25]]) @ transition_19_3









    Out[20]:





array([[0.358, 0.642]])

Exercise 1.19 c)



In [21]:

    
stationary_distribution(transition_19)









    Out[21]:





array([0.25, 0.75])

Exercise 1.20



In [22]:

    
transition_20 = np.array([[0.75, 0.25],
                          [0.08, 0.92]])
v_20_1995 = np.array([[0.7, 0.3]])

Exercise 1.20 a)



In [23]:

    
transition_20_3 = np.linalg.matrix_power(transition_20, 3)
transition_20_3









    Out[23]:





array([[0.470275, 0.529725],
       [0.169512, 0.830488]])



In [24]:

    
v_20_1995 @ transition_20_3









    Out[24]:





array([[0.3800461, 0.6199539]])

Exercise 1.20 b)



In [25]:

    
transition_20_10 = np.linalg.matrix_power(transition_20, 10)
transition_20_10









    Out[25]:





array([[0.25623362, 0.74376638],
       [0.23800524, 0.76199476]])



In [26]:

    
v_20_1995 @ transition_20_10









    Out[26]:





array([[0.25076511, 0.74923489]])

Exercise 1.20 c)



In [27]:

    
stationary_distribution(transition_20)









    Out[27]:





array([0.24242424, 0.75757576])

Exercise 1.36



In [28]:

    
transition_36 = np.array([[   0,    0,    1],
                          [0.05, 0.95,    0],
                          [   0, 0.02, 0.98]])

Exercise 1.36 a)



In [29]:

    
stationary_distribution(transition_36)









    Out[29]:





array([0.01408451, 0.28169014, 0.70422535])

Exercise 1.36 b)



In [30]:

    
1 / stationary_distribution(transition_36)









    Out[30]:





array([71.  ,  3.55,  1.42])

Exercise 1.37



In [31]:

    
transition_37 = np.array([[  1,   0,   0,   0],
                          [0.2, 0.8,   0,   0],
                          [  0, 0.2, 0.8,   0],
                          [  0,   0, 0.2, 0.2]])

Exercise 1.37 a)



In [32]:

    
stationary_distribution(transition_37)









    Out[32]:





array([1., 0., 0., 0.])

Exercise 1.38



In [33]:

    
transition_38 = np.array([[0.5, 0.5,   0,   0],
                          [2/3,   0, 1/3,   0],
                          [3/4,   0,   0, 1/4],
                          [  1,   0,   0,   0]])

Exercise 1.38 a)



In [34]:

    
pi_38 = stationary_distribution(transition_38)

Exercise 1.38 b)



In [35]:

    
check_detailed_balance(transition_38)









    Out[35]:





False

Exercise 1.48



In [36]:

    
n = 12
transition_48 = np.zeros([n, n])
for i in range(n):
    transition_48[i, [(i - 1) % n, (i + 1) % n]] = [0.5] * 2
transition_48









    Out[36]:





array([[0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5],
       [0.5, 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0.5, 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 0.5, 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0.5, 0. , 0.5, 0. , 0. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0. , 0.5, 0. , 0.5, 0. , 0. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 0.5, 0. , 0.5, 0. , 0. , 0. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0. , 0.5, 0. , 0. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0. , 0.5, 0. , 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0. , 0.5, 0. ],
       [0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0. , 0.5],
       [0.5, 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0. , 0.5, 0. ]])

Exercise 1.48 a)

El resultado será: $$\pi_{i} = 1 / 12 \ \forall i \in {1, 2, ..., 12}$$ Por cumplir la matriz de transición la propiedad de ser doblemente estocástica (tanto filas como columnas suman la unidad).



In [37]:

    
stationary_distribution(transition_48)









    Out[37]:





array([0.08333333, 0.08333333, 0.08333333, 0.08333333, 0.08333333,
       0.08333333, 0.08333333, 0.08333333, 0.08333333, 0.08333333,
       0.08333333, 0.08333333])

Exercise 1.48 b)



In [39]:

    
import numpy as np

def position_to_clock(position: int) -> int:
    if (position % 12):
        return position % 12
    else:
        return 12

n = 100000
y = 0
for n_temp in range(1, n + 1):
    s = f'{12} '
    position = np.random.choice([-1, 1])
    s += f'{position_to_clock(position): 2d} '
    while(position % 12 != 0):
        position += np.random.choice([-1, 1])
        s += f'{position_to_clock(position): 2d} '
    y += (position != 0)
    # print(f'y: {y:3d} s: {s}')
    if n_temp % 10000 == 0:
        print(f'n = {n_temp}\tp = {y / n_temp:.5f}')
print(f'n = {n}\tp = {y / n:.5f}')









    



n = 10000	p = 0.08410
n = 20000	p = 0.08265
n = 30000	p = 0.08237
n = 40000	p = 0.08260
n = 50000	p = 0.08374
n = 60000	p = 0.08283
n = 70000	p = 0.08313
n = 80000	p = 0.08273
n = 90000	p = 0.08238
n = 100000	p = 0.08262
n = 100000	p = 0.08262