import numpy as np
import matplotlib.pyplot as plt

N = 10000 # nombre de tirages
X = ??? # simuler un vecteur de v.a. uniformes sur [0,1]

fig = plt.figure()
plt.hist(X, density="True", bins=100, label="histogramme des realisations")
plt.plot([0, 1], [1, 1], color="red")
plt.show()

## Simuler une loi exponentielle en utilisant np.random.exponential

N = 10000 # nombre de tirages
lam = 2 # parametre de la loi exponentielle
Z = ??? # utiliser la fonction np.random.exponential

x = np.linspace(0,10,100)
densiteExp = lam * np.exp(-lam*x)

fig = plt.figure()
plt.hist(Z, density="True", bins=100, label="histogramme des realisations")
plt.plot(x, densiteExp, color="red", label="densite")
plt.show()

## Simuler une loi exponentielle a partir d'une loi uniforme

N = 10000 # nombre de tirages
Z = ??? # utiliser un vecteur de v.a. uniformes sur [0,1] et l'inversion de la fonction de repartition

fig = plt.figure()
plt.hist(Z, density="True", bins=100, label="histogramme des realisations")
plt.plot(x, densiteExp, color="red", label="densite")
plt.show()

## Simuler deux variables aleatoires gaussiennes independantes avec la methode de Box-Muller

N = 10000 # nombre de paires de variables gaussiennes
U1 = np.random.rand(N) # premiere variable uniforme sur [0,1]
U2 = np.random.rand(N) # deuxieme variable uniforme sur [0,1]

# Methode de Box-Muller : si U1 et U2 sont independantes et uniformes sur (0,1],
# alors Z1 = sqrt(-2*log(U1))*cos(2*pi*U2) et Z2 = sqrt(-2*log(U1))*sin(2*pi*U2)
# sont deux variables aleatoires gaussiennes N(0,1) independantes

Z1 = ???
Z2 = ???

x = np.linspace(-4, 4, 100)
densiteGauss = (1 / np.sqrt(2 * np.pi)) * np.exp(-x**2 / 2)

fig, axes = plt.subplots(1, 2, figsize=(12, 5))

# Histogramme de Z1
axes[0].hist(Z1, density=True, bins=100, label="histogramme de Z1")
axes[0].plot(x, densiteGauss, color="red", label="densite N(0,1)")
axes[0].set_title("Distribution de Z1 (Box-Muller)")
axes[0].legend()

# Histogramme de Z2
axes[1].hist(Z2, density=True, bins=100, label="histogramme de Z2")
axes[1].plot(x, densiteGauss, color="red", label="densite N(0,1)")
axes[1].set_title("Distribution de Z2 (Box-Muller)")
axes[1].legend()

plt.tight_layout()
plt.show()

# Verifier l'independance : correlation entre Z1 et Z2 devrait etre proche de 0
correlation = np.corrcoef(Z1, Z2)[0, 1]
print(f"Correlation entre Z1 et Z2: {correlation:.6f}")

# Scatter plot de (Z1, Z2) pour visualiser l'independance
fig = plt.figure(figsize=(8, 8))
plt.scatter(Z1, Z2, alpha=0.3, s=10)
plt.xlabel('Z1')
plt.ylabel('Z2')
plt.title('Scatter plot de (Z1, Z2) - Variables gaussiennes independantes (Box-Muller)')
plt.grid(True, alpha=0.3)
plt.axis('equal')
plt.show()

## Generer un vecteur gaussien general en dimension 2 a partir de (Z1, Z2)

# Parametres du vecteur gaussien general N(mu, Gamma)
mu = np.array([1, 2])  # moyenne
sigma1 = 0.5  # ecart-type de la premiere composante
sigma2 = 1.5  # ecart-type de la deuxieme composante
rho = 0.7     # correlation entre les deux composantes

# Matrice de covariance
# Gamma = [[sigma1^2, rho*sigma1*sigma2], [rho*sigma1*sigma2, sigma2^2]]
cov_matrix = np.array([
    [sigma1**2, rho * sigma1 * sigma2],
    [rho * sigma1 * sigma2, sigma2**2]
])

# Decomposition de Cholesky: Gamma = L * L^T
L = np.linalg.cholesky(cov_matrix)

# Transformer (Z1, Z2) en un vecteur gaussien general
# X = mu + L * [Z1, Z2]^T
Z = np.column_stack([Z1, Z2])  # matrice N x 2
X = mu + Z @ L.T  # transformation affine

X1 = X[:, 0]
X2 = X[:, 1]

# Affichage des resultats
fig, axes = plt.subplots(2, 2, figsize=(12, 10))

# Histogramme de X1 - top left
axes[0, 0].hist(X1, density=True, bins=100, label="histogramme de X1")
x1_range = np.linspace(X1.min(), X1.max(), 200)
densiteX1 = (1 / (sigma1 * np.sqrt(2 * np.pi))) * np.exp(-(x1_range - mu[0])**2 / (2 * sigma1**2))
axes[0, 0].plot(x1_range, densiteX1, color="red", linewidth=2, label="densite N(mu1, sigma1)")
axes[0, 0].set_title(f"Distribution de X1 (mu={mu[0]}, sigma={sigma1})")
axes[0, 0].set_xlabel('X1')
axes[0, 0].legend()

# Matrice de correlation empirique - top right
correlation_empirique = np.corrcoef(X1, X2)[0, 1]
axes[0, 1].text(0.1, 0.8, f'Parametres du vecteur gaussien:', fontsize=12, fontweight='bold')
axes[0, 1].text(0.1, 0.7, f'mu = {mu}', fontsize=11)
axes[0, 1].text(0.1, 0.5, f'Matrice de covariance:\n{cov_matrix}', fontsize=10)
axes[0, 1].text(0.1, 0.3, f'Correlation theorique: {rho:.3f}', fontsize=11)
axes[0, 1].text(0.1, 0.2, f'Correlation empirique: {correlation_empirique:.3f}', fontsize=11)
axes[0, 1].axis('off')

# Scatter plot de (X1, X2) - bottom left
axes[1, 0].scatter(X1, X2, alpha=0.3, s=10)
axes[1, 0].set_xlabel('X1')
axes[1, 0].set_ylabel('X2')
axes[1, 0].set_title(f'Scatter plot de (X1, X2) - Vecteur gaussien avec rho={rho}')
axes[1, 0].grid(True, alpha=0.3)

# Histogramme de X2 (rotated 90 degrees to the right) - bottom right
axes[1, 1].hist(X2, density=True, bins=100, label="histogramme de X2", orientation="horizontal")
x2_range = np.linspace(X2.min(), X2.max(), 200)
densiteX2 = (1 / (sigma2 * np.sqrt(2 * np.pi))) * np.exp(-(x2_range - mu[1])**2 / (2 * sigma2**2))
axes[1, 1].plot(densiteX2, x2_range, color="red", linewidth=2, label="densite N(mu2, sigma2)")
axes[1, 1].set_title(f"Distribution de X2 (mu={mu[1]}, sigma={sigma2})")
axes[1, 1].set_ylabel('X2')
axes[1, 1].legend()

plt.tight_layout()
plt.show()

print(f"\nVecteur gaussien general N(mu, Gamma):")
print(f"Moyenne empirique: [{X1.mean():.4f}, {X2.mean():.4f}]")
print(f"Ecart-types empiriques: [{X1.std():.4f}, {X2.std():.4f}]")
print(f"Correlation empirique: {correlation_empirique:.4f}")