Simulation code to solve the Riddler Classic puzzle (Supreme Court Judges),https://fivethirtyeight.com/features/how-many-bingo-cards-are-there-in-the-world/

main.py

import numpy as np

def runSimulation(N):
    judge_num = np.zeros(N, dtype=np.int)

    for i in range(9):
        judge_num[:(np.random.randint(0, 41))] += 1

    for i in range(4, N, 4):
        if np.random.random() > 0.5:
            # the same party has the presidency and congress so begin filling up the seat
            while judge_num[i] < 9:
                judge_num[i:i + (np.random.randint(0, 41))] += 1

    print(9 - np.mean(judge_num))


# run main simulation
runSimulation(400*1000*1000)
runSimulation(400*1000*1000)
runSimulation(400*1000*1000)

Expected number of vacancies is 1.904