chaincode.py

import numpy as np

def HHH(p):
  return 1/p + 1/(p**2) + 1/(p**3)

tosses = []
for i in range(10000):
  chain = []
  while True:
    chain.append(np.random.choice([0, 1], p=[.5, .5]))
    if len(chain) > 2:
      if chain[-1] == 1 and chain[-2] == 1 and chain[-3] == 1:
        tosses.append(len(chain))
        break
print(HHH(1/2))
print(np.mean(tosses))