

In [ ]:



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.pinv(A)[len(P) - 1, :]
    return p_stationary

Exercise: Ruiz Family

La familia Ruiz recibe el periódico todas las mañanas, y lo coloca en el revistero después de leerlo. Cada tarde, con probabilidad 0.3, alguien coge todos los periódicos del revistero y los lleva al contenedor de papel. Por otro lado, si hay al menos 5 periódicos en el montón, el señor Ruiz los lleva al contenedor.

a) Construye una cadena de Markov que cuente el número de periódicos que hay en el revistero cada noche. ¿Cómo son los estados?



In [3]:

    
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]])

b) Si el domingo por la noche está vacío el revistero, ¿Cuál es la probabilidad de que haya 1 periódico el miércoles por la noche?



In [4]:

    
np.linalg.matrix_power(transition_ruiz, 4)









    Out[4]:





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.   ]])

c) Calcula la probabilidad, a largo plazo, de que el revistero esté vacío una noche cualquiera.



In [5]:

    
stationary_distribution(transition_ruiz)









    Out[5]:





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

Exercise 1.36



In [6]:

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

Exercise 1.36 a)



In [7]:

    
stationary_distribution(transition_36)









    Out[7]:





array([0.01408451, 0.28169014, 0.70422535])

Exercise 1.36 b)



In [8]:

    
1 / stationary_distribution(transition_36)









    Out[8]:





array([71.  ,  3.55,  1.42])

Exercise 1.48



In [9]:

    
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[9]:





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)

La distribución estacionaria 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). Por tanto, dado que se pide el número medio de pasos: $$ E_i(T_i) = \frac{1}{\pi_i} = \frac{1}{1 / 12} = 12 \ \forall i \in {1, 2, ..., 12}$$

El mismo resultado se obtiene al ejecutar las operaciones:



In [10]:

    
1 / stationary_distribution(transition_48)









    Out[10]:





array([12., 12., 12., 12., 12., 12., 12., 12., 12., 12., 12., 12.])

Exercise 1.48 b)



In [11]:

    
import numpy as np

n = 100000
y = 0
d = 12
for n_temp in range(1, n + 1):
    visited = set()
    k = np.random.choice(range(d))
    
    position = k
    s = str(position % d) + ' '
    visited.add(position % d)
    
    position += np.random.choice([-1, 1])
    s += str(position % d) + ' '
    visited.add(position % d)
    
    while(position % d != k):
        position += np.random.choice([-1, 1])
        s += str(position % d) + ' '
        visited.add(position % d)
        
    y += (len(visited) == d)
    
    # print(y, s, sep=', ')
    if n_temp % 10000 == 0:
        print(y / n_temp)
print(y / n)









    



0.0881
0.08825
0.09096666666666667
0.0914
0.09116
0.09143333333333334
0.0909
0.0906625
0.09092222222222222
0.09117
0.09117