import random
import numpy as np
import matplotlib.pyplot as plt

# Paramètres
b = 2 #nb de boules dans l'urne de départ
N = 2000 # nb de simulations
T = 30 #nb d'échanges

# X[k, n] : nombre de boules dans l'urne 1
X = np.zeros((N, T + 1), dtype=int)

for k in range(N):
    # 0 : urne 2 ; 1 : urne 1
    urne = [0] * b

    for n in range(T):
        j = random.randrange(b)
        urne[j] = 1 - urne[j]
        X[k, n + 1] = sum(urne)

# Moyenne expérimentale
moyenne = X.mean(axis=0)

# Espérance théorique
n = np.arange(T + 1)
theorie = b / 2 * (1 - (1 - 2 / b) ** n)

# Graphiques
fig, axes = plt.subplots(1, 2, figsize=(13, 5))

# Une trajectoire
axes[0].plot(n, X[0], linewidth=1)
axes[0].axhline(b / 2, linestyle="--",
                label="b/2")
axes[0].set_title("Une simulation")
axes[0].set_xlabel("Instant n")
axes[0].set_ylabel("Nombre de boules dans l'urne 1")
axes[0].legend()
axes[0].grid()

# Comparaison des espérances
axes[1].plot(n, moyenne, label="Moyenne simulée")
axes[1].plot(n, theorie, "--", label="Espérance théorique")
axes[1].set_title(f"Moyenne sur {N} simulations")
axes[1].set_xlabel("Instant n")
axes[1].set_ylabel("E(X_n)")
axes[1].legend()
axes[1].grid()

plt.tight_layout()
plt.show()