# %% [markdown]
# Validation of SBC Model with Observational $H(z)$ Data

This notebook reproduces the validation of the SBC model using observational $H(z)$ data. We perform both direct parameter fitting and MCMC posterior sampling to estimate the best-fit values of $H_0$ and $q_0$, followed by residual analysis and statistical metrics.


# %%
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import emcee
from scipy.optimize import curve_fit
from scipy.stats import chi2
import corner
import os

plt.style.use("ggplot")
%matplotlib inline

# %% [markdown]
## Load Observational $H(z)$ Data

# %%
file_path = "../data/observational_Hz.csv"
if os.path.exists(file_path):
    df = pd.read_csv(file_path)
    z_data = df["z"].values
    Hz_data = df["Hz"].values
    Hz_err = df["Hz_err"].values
    print(f"✅ Loaded {len(df)} H(z) data points.")
else:
    raise FileNotFoundError(f"❌ File not found at {file_path}. Please check the relative path.")

# %% [markdown]
## Define SBC Expansion Model

# %%
c = 299792.458  # speed of light in km/s

def Hz_SBC(z, H0, q0):
    return H0 * (1 + q0 * z)

# %% [markdown]
## Least Squares Fit

# %%
popt_int, pcov_int = curve_fit(Hz_SBC, z_data, Hz_data, sigma=Hz_err, absolute_sigma=True)
H0_fit, q0_fit = popt_int
H0_err, q0_err = np.sqrt(np.diag(pcov_int))

print(f"H0 = {H0_fit:.2f} ± {H0_err:.2f} km/s/Mpc")
print(f"q0 = {q0_fit:.2f} ± {q0_err:.2f}")

# %% [markdown]
## Plot Fit and Residuals

# %%
z_plot = np.linspace(0, 2, 500)
Hz_model = Hz_SBC(z_plot, H0_fit, q0_fit)

fig, ax = plt.subplots(figsize=(8, 5))
ax.errorbar(z_data, Hz_data, yerr=Hz_err, fmt='o', label="Observed H(z)", color="black")
ax.plot(z_plot, Hz_model, label=f"SBC Fit: $H_0$ = {H0_fit:.2f}, $q_0$ = {q0_fit:.2f}", color="red")
ax.set_xlabel("Redshift $z$")
ax.set_ylabel("$H(z)$ [km/s/Mpc]")
ax.legend()
ax.grid(True)
plt.show()
