import numpy as np
import matplotlib.pyplot as plt

def bernoulli_density(x, p = 0.5):
    return p ** x * (1 - p) ** (1 - x)

bernoulli_density(0, p = 0.27)

0.73

bernoulli_density(1, p = 0.27)

0.27

rng = np.random.default_rng()
x = rng.binomial(1, 0.27, size = 10)
bernoulli_density(x, p = 0.27)

array([0.73, 0.73, 0.73, 0.73, 0.73, 0.73, 0.73, 0.27, 0.27, 0.73])

x = np.array([0, 1])
p = 0.33
y = bernoulli_density(x, p = p)
plt.scatter(x, y)

N = 100000
R = 10100
x = rng.binomial(1, p, size = (R, N))
ndx = np.arange(1, N + 1)
cm = np.cumsum(x, axis = 1) / ndx

plt.plot(cm.T);
plt.axhline(p, color = "black", linestyle = "--");

np.shape(cm)

(10100, 100000)

plt.hist(cm[:, 800]);