# 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