"""
Module: Regression Lineaire Bayesienne
Categorie: Bayesian Statistics
Difficulte: Intermediaire

Genere depuis la plateforme ML Formation
"""

# Imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from sklearn.model_selection import train_test_split
from sklearn.metrics import accuracy_score, mean_squared_error, r2_score

# Charger le dataset
df = pd.read_csv('housing_simple.csv')

# Rappel OLS et preparation des donnees
# Type: Code executable
from sklearn.linear_model import LinearRegression

print("=" * 70)
print("       POINT DE DEPART: LA REGRESSION CLASSIQUE (OLS)")
print("=" * 70)

print("""
On reprend le dataset immobilier du module Regression Lineaire:
predire le prix a partir de la surface. Pour que les priors gaussiens
centres sur 0 aient du sens, on STANDARDISE les variables (moyenne 0,
ecart-type 1): une pratique quasi obligatoire en bayesien.
""")

display(df.head(8), title="Dataset immobilier")

# Variables standardisees (z-scores)
x_raw = df["surface"].values.astype(float)
y_raw = df["price"].values.astype(float)
x = (x_raw - x_raw.mean()) / x_raw.std()
y = (y_raw - y_raw.mean()) / y_raw.std()
n = len(x)

# OLS de reference avec sklearn
ols = LinearRegression()
ols.fit(x.reshape(-1, 1), y)
a_ols = float(ols.coef_[0])
b_ols = float(ols.intercept_)

# Estimation du bruit d'observation sigma2 (variance des residus OLS)
residus = y - ols.predict(x.reshape(-1, 1))
sigma2 = float(np.var(residus, ddof=2))

print(f"Nombre de ventes        : {n}")
print(f"Pente OLS (standardisee): {a_ols:.4f}")
print(f"Intercept OLS           : {b_ols:.4f}")
print(f"Variance du bruit sigma2: {sigma2:.4f}")

print("\n" + "-" * 40)
print("LECTURE:")
print("-" * 40)
print(f"""
• En unites standardisees, la pente {a_ols:.3f} signifie: +1 ecart-type
  de surface -> +{a_ols:.3f} ecart-type de prix. La correlation est forte.
• sigma2 = {sigma2:.3f} est la part de variance NON expliquee par la droite:
  c'est le bruit que le modele bayesien devra integrer.
• L'OLS donne ces nombres... et rien d'autre. Combien vaut l'incertitude
  sur la pente ? Motus. La suite du module la calcule exactement.
""")

# Nuage de points avec la droite OLS
fig, ax = plt.subplots(figsize=(9, 5.5))
ax.scatter(x, y, color="#9B7AC4", alpha=0.7, label="Ventes (standardisees)")
grille = np.linspace(x.min(), x.max(), 100)
ax.plot(grille, a_ols * grille + b_ols, color="#e67e22", lw=2.5,
        label=f"Droite OLS (pente {a_ols:.3f})")
ax.set_xlabel("Surface (z-score)")
ax.set_ylabel("Prix (z-score)")
ax.set_title("La regression classique: une droite unique, zero incertitude affichee")
ax.legend()
plt.tight_layout()
plt.show()


# Le posterior exact des coefficients
# Type: Code executable
print("=" * 70)
print("       POSTERIOR ANALYTIQUE DES COEFFICIENTS (a, b)")
print("=" * 70)

print("""
Modele: y = a*x + b + bruit, bruit ~ N(0, sigma2)
Prior : (a, b) ~ N(0, tau2 * I)

La conjugaison gaussienne donne un posterior gaussien EXACT:
  Cov_post = inverse( X'X / sigma2 + I / tau2 )
  mu_post  = Cov_post @ X'y / sigma2

(X est la matrice [x, 1] pour porter la pente ET l'intercept.)
""")

# Preparation (identique a la cellule precedente)
x_raw = df["surface"].values.astype(float)
y_raw = df["price"].values.astype(float)
x = (x_raw - x_raw.mean()) / x_raw.std()
y = (y_raw - y_raw.mean()) / y_raw.std()
n = len(x)
sigma2 = 0.05  # bruit d'observation (proche de l'estimation OLS)
tau2 = 1.0     # prior N(0, 1) sur chaque coefficient

# Matrice de design [x, 1]
X = np.column_stack([x, np.ones(n)])

# Formules du posterior gaussien
precision_post = X.T @ X / sigma2 + np.eye(2) / tau2
cov_post = np.linalg.inv(precision_post)
mu_post = cov_post @ X.T @ y / sigma2

a_mean, b_mean = mu_post
a_std = float(np.sqrt(cov_post[0, 0]))
b_std = float(np.sqrt(cov_post[1, 1]))

print(f"Posterior de la pente a    : moyenne {a_mean:.4f}, ecart-type {a_std:.4f}")
print(f"Posterior de l'intercept b : moyenne {b_mean:.4f}, ecart-type {b_std:.4f}")
print(f"Intervalle de credibilite 95% pour a : "
      f"[{a_mean - 1.96 * a_std:.4f}, {a_mean + 1.96 * a_std:.4f}]")

# Densites marginales des coefficients
fig, axes = plt.subplots(1, 2, figsize=(11, 4.5))
for ax, m, s, nom, couleur in [
    (axes[0], a_mean, a_std, "pente a", "#27ae60"),
    (axes[1], b_mean, b_std, "intercept b", "#3498db"),
]:
    grille = np.linspace(m - 4 * s, m + 4 * s, 400)
    ax.plot(grille, stats.norm(m, s).pdf(grille), color=couleur, lw=2.5)
    ax.fill_between(grille, stats.norm(m, s).pdf(grille), alpha=0.15, color=couleur)
    ax.axvline(m, color=couleur, linestyle="--")
    ax.set_title(f"Posterior de la {nom}")
    ax.set_xlabel(nom)
    ax.set_ylabel("Densite")
plt.tight_layout()
plt.show()

print("\n" + "-" * 40)
print("INTERPRETATION:")
print("-" * 40)
print(f"""
• La pente n'est plus un nombre mais une GAUSSIENNE: centree sur
  {a_mean:.3f}, avec une vraie barre d'erreur ({a_std:.3f}).
• Le prior N(0, 1) tire legerement les coefficients vers 0 par rapport
  a l'OLS: c'est l'effet regularisant (Ridge!) du prior.
• Essayez tau2 = 0.01 (prior tres serre) ou sigma2 plus grand (donnees
  plus bruitees) et observez comme le posterior se deplace: le bras de
  fer prior/donnees de l'exercice manuel, en vrai.
""")


# Visualiser l'incertitude: le faisceau de droites plausibles
# Type: Code executable
print("=" * 70)
print("       50 DROITES TIREES DU POSTERIOR")
print("=" * 70)

print("""
Le posterior des coefficients est une distribution sur (a, b)...
donc une distribution sur des DROITES ENTIERES. En echantillonner 50,
c'est visualiser directement "toutes les regressions compatibles
avec les donnees et le prior".
""")

# Preparation et posterior (identiques a la cellule precedente)
x_raw = df["surface"].values.astype(float)
y_raw = df["price"].values.astype(float)
x = (x_raw - x_raw.mean()) / x_raw.std()
y = (y_raw - y_raw.mean()) / y_raw.std()
n = len(x)
sigma2, tau2 = 0.05, 1.0
X = np.column_stack([x, np.ones(n)])
cov_post = np.linalg.inv(X.T @ X / sigma2 + np.eye(2) / tau2)
mu_post = cov_post @ X.T @ y / sigma2

# Echantillonner 50 couples (a, b) du posterior gaussien bivarie
rng = np.random.default_rng(42)
tirages = rng.multivariate_normal(mu_post, cov_post, size=50)

grille = np.linspace(x.min() - 0.3, x.max() + 0.3, 100)
fig, ax = plt.subplots(figsize=(10, 6))
for a_i, b_i in tirages:
    ax.plot(grille, a_i * grille + b_i, color="#C09CF0", alpha=0.25, lw=1)
ax.plot(grille, mu_post[0] * grille + mu_post[1], color="#9B7AC4", lw=2.5,
        label="Droite moyenne du posterior")
ax.scatter(x, y, color="#3A3A3A", zorder=5, alpha=0.8, label="Ventes observees")
ax.set_xlabel("Surface (z-score)")
ax.set_ylabel("Prix (z-score)")
ax.set_title("Le faisceau des regressions plausibles (50 tirages du posterior)")
ax.legend()
plt.tight_layout()
plt.show()

# Ou le faisceau s'ecarte-t-il le plus ?
ecarts = np.array([[a_i * g + b_i for g in grille] for a_i, b_i in tirages]).std(axis=0)
print(f"Dispersion du faisceau au centre du nuage : {ecarts[len(grille) // 2]:.4f}")
print(f"Dispersion du faisceau aux extremites     : {ecarts[0]:.4f} / {ecarts[-1]:.4f}")

print("\n" + "-" * 40)
print("CE QU'IL FAUT VOIR:")
print("-" * 40)
print("""
• Toutes les droites se pincent au CENTRE du nuage (la, les donnees
  contraignent fort) et s'evasent aux EXTREMITES (la, on extrapole).
• L'incertitude n'est pas un nombre unique: elle depend de l'endroit
  ou l'on predit. La cellule suivante la transforme en bandes de
  prediction exploitables.
""")


# Predictions avec incertitude: les bandes de credibilite
# Type: Code executable
from sklearn.linear_model import LinearRegression

print("=" * 70)
print("       BANDES DE CREDIBILITE PREDICTIVES A 95%")
print("=" * 70)

print("""
Pour un point x*, la prediction bayesienne est une GAUSSIENNE:
  moyenne  : mu_post . x*
  variance : x*' Cov_post x*   (incertitude sur les coefficients)
             + sigma2           (bruit d'observation irreductible)

On compare avec la prediction ponctuelle de sklearn.
""")

# Preparation et posterior (identiques aux cellules precedentes)
x_raw = df["surface"].values.astype(float)
y_raw = df["price"].values.astype(float)
x = (x_raw - x_raw.mean()) / x_raw.std()
y = (y_raw - y_raw.mean()) / y_raw.std()
n = len(x)
sigma2, tau2 = 0.05, 1.0
X = np.column_stack([x, np.ones(n)])
cov_post = np.linalg.inv(X.T @ X / sigma2 + np.eye(2) / tau2)
mu_post = cov_post @ X.T @ y / sigma2

# Prediction bayesienne sur une grille
grille = np.linspace(x.min() - 0.5, x.max() + 0.5, 120)
X_star = np.column_stack([grille, np.ones(len(grille))])
pred_mean = X_star @ mu_post
var_coef = np.sum(X_star @ cov_post * X_star, axis=1)  # x*' Cov x* par ligne
pred_std = np.sqrt(var_coef + sigma2)

# Prediction ponctuelle sklearn (OLS)
ols = LinearRegression().fit(x.reshape(-1, 1), y)
pred_sklearn = ols.predict(grille.reshape(-1, 1))

fig, ax = plt.subplots(figsize=(10, 6))
ax.fill_between(grille, pred_mean - 1.96 * pred_std, pred_mean + 1.96 * pred_std,
                color="#C09CF0", alpha=0.30, label="Bande de credibilite 95%")
ax.plot(grille, pred_mean, color="#9B7AC4", lw=2.5, label="Prediction bayesienne (moyenne)")
ax.plot(grille, pred_sklearn, color="#e67e22", lw=1.8, linestyle="--",
        label="Prediction sklearn (ponctuelle)")
ax.scatter(x, y, color="#3A3A3A", alpha=0.8, zorder=5, label="Ventes observees")
ax.set_xlabel("Surface (z-score)")
ax.set_ylabel("Prix (z-score)")
ax.set_title("La meme droite... plus l'incertitude qui va avec")
ax.legend()
plt.tight_layout()
plt.show()

# Exemple chiffre au centre et en extrapolation
for x_star in [0.0, grille.max()]:
    v = np.array([x_star, 1.0])
    m = float(v @ mu_post)
    s = float(np.sqrt(v @ cov_post @ v + sigma2))
    zone = "centre du nuage" if x_star == 0.0 else "extrapolation  "
    print(f"  x* = {x_star:+.2f} ({zone}) : prix predit {m:+.3f} "
          f"± {1.96 * s:.3f} (95%)")

print("\n" + "-" * 40)
print("INTERPRETATION:")
print("-" * 40)
print("""
• Les deux courbes moyennes sont presque confondues: sur ce dataset,
  bayesien et OLS s'accordent sur la tendance.
• La VRAIE difference est la bande violette: plus etroite au centre,
  plus large en extrapolation, exactement comme le faisceau de droites.
• En pratique (pricing, medecine, industrie), cette bande change tout:
  elle dit QUAND le modele sait, et quand il devine.
""")


# Exercice: le prior contre Ridge, match retour
# Type: Exercice
# Exercice: verifier numeriquement que MAP bayesien = Ridge
# avec la correspondance alpha = sigma2 / tau2

from sklearn.linear_model import Ridge

# Preparation standard
x_raw = df["surface"].values.astype(float)
y_raw = df["price"].values.astype(float)
x = (x_raw - x_raw.mean()) / x_raw.std()
y = (y_raw - y_raw.mean()) / y_raw.std()
n = len(x)
X = np.column_stack([x, np.ones(n)])
sigma2 = 0.05

# TODO: pour tau2 dans [10.0, 1.0, 0.1, 0.01]:
#   1. calculez la moyenne du posterior bayesien mu_post
#      (formules de la cellule "posterior exact")
#   2. calculez alpha = sigma2 / tau2 et ajustez
#      Ridge(alpha=alpha, fit_intercept=False) sur (X, y)
#      (fit_intercept=False car l'intercept est deja dans X)
#   3. comparez la pente bayesienne et la pente Ridge
# TODO: que devient la pente quand tau2 retrecit ? Pourquoi ?