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