Bernoulli two pebbles patch.py

From Werner KRAUTH

import random
p = 0.7
phat = 0.4
count_a = 0
N_trial = 1000000
for iter in range(N_trial):
    Upsilon_1 = random.uniform(0.0, 1.0)
    if Upsilon_1 < phat: count_a += 1
    else:
        Upsilon_2 = random.uniform(0.0, 1.0)
        if Upsilon_2 > (1 - p) / (1 - phat): count_a += 1
print(count_a / N_trial)