markovcoinflip.py

import numpy as np

N = 10

P0 = np.array([[.5, .5, 0, 0], [.5, 0, .5, .0], 
[.5, 0, 0, .5], [.5, 0, 0, .5]])

P1 = np.array([[.5, .5, 0, 0], [.5, 0, .5, .0], 
[.5, 0, 0, .5], [0, 0, 0, 1]])

hits = []

# Simulation 1
for i in range(10000):
    chain = []
    flag = True
    for j in range(N):
        X = np.random.choice([0, 1], p=[1/2, 1/2])
        chain.append(X)
        if len(chain)>2:
            if chain[-1] == 1 and chain[-2] == 1 and chain[-3] == 1:
                hits.append(1)
                flag = False
                break
    if flag:
        hits.append(0)

hits1 = []
# Simulation 2
for i in range(10000):
    chain = []
    flag = True
    for j in range(N):
        X = np.random.choice([0, 1], p=[1/2, 1/2])
        chain.append(X)
    if chain[-1] == 1 and chain[-2] == 1 and chain[-3] == 1:
        hits1.append(1)
    else:
      hits1.append(0)

hits2 = []
# Simulation 3
for i in range(10000):
    chain = []
    chains = 0
    flag = True
    for j in range(N):
        X = np.random.choice([0, 1], p=[1/2, 1/2])
        chain.append(X)
        if len(chain)>2:
            if chain[-1] == 1 and chain[-2] == 1 and chain[-3] == 1:
                chains += 1
    if chain[-1] == 1 and chain[-2] == 1 and chain[-3] == 1 and chains == 1:
        hits2.append(1)
    else:
      hits2.append(0)

print('Simulated Probability of Seeing 3H Chain in N Tosses', np.mean(hits))
print('Simulated Probability of Seeing 3H Chain at the End of N Tosses', np.mean(hits1))
print('Simulated Probability of Seeing 3H Chain for the first time at the End of N Tosses', np.mean(hits2))
print('Q1', np.linalg.matrix_power(P1, N)[0][3])
print('Q2', np.linalg.matrix_power(P0, N)[0][3])
print('Q3', np.linalg.matrix_power(P1, N)[0][3] - np.linalg.matrix_power(P1, N-1)[0][3])