coin_flip_sim.py

import numpy as np

# transition matrix
P = np.array([[.5, .5, 0, 0],
              [.5, 0, .5, 0],
              [.5, 0, 0, .5],
              [.5, 0, 0, .5]])

# Simulation
hits = []
for i in range(10000):
  chain = []
  for j in range(10):
    chain.append(np.random.choice([0, 1], p=[.5, .5]))
  if chain[-1] == 1 and chain[-2] == 1 and chain[-3] == 1:
    hits.append(1)
  else:
    hits.append(0)
  

print('Simulated Probability: ', np.mean(hits))
print('Analytical Probability: ', np.linalg.matrix_power(P, 10)[0][3])