ckprb.py

import numpy as np


P = np.array([[.7, .3],[.6, .4]])

P26 = np.linalg.matrix_power(P, 26)

print(P26[1][0]*.5 + P26[0][0]*.5)

day27 = []
for i in range(10000):
    
    # Randomly select coin 1 or 2 on day 1
    X = np.random.randint(0, 2)
    
    for j in range(26):
        # Flip a coin for 26 days
        if X == 1:
            X = np.random.choice([0, 1], p=[.3, .7])
        else:
            X = np.random.choice([0, 1], p=[.4, .6])
    
    if X == 1:
        day27.append(1)
    else:
        day27.append(0)

print(np.mean(day27))