04 LPU Ensemble – CI Band Visualization#
Local Parameter Uncertainty via SciPy leastsq Jacobian.
Fits Q_sim = a*WAM + b, extracts parameter covariance from inv(JtJ),
and propagates uncertainty via MVN sampling of [a, b] parameter pairs.
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from pathlib import Path
from red_tide_reanalysis.ingestion.obs_loader import load_observations
from red_tide_reanalysis.ingestion.wam_loader import load_wam_model
from red_tide_reanalysis.ingestion.align import align_obs_model
from red_tide_reanalysis.lpu.method import LPUMethod
from red_tide_reanalysis.writers.ensemble_writer import write_ensemble_csv
from red_tide_reanalysis.writers.stats_writer import write_stats_csv
1. Load and Align Data#
obs_raw = load_observations("../Observation_Data/Observed_flow_ARCADIA_FL.csv")
wam_raw = load_wam_model("../Synthetic_Model_Data/Station 02296750 (ARCADIA)_reach000084_83.csv")
obs, model = align_obs_model(obs_raw, wam_raw)
print(f"Aligned series: {len(obs)} timesteps, {obs.index[0].date()} to {obs.index[-1].date()}")
Aligned series: 9131 timesteps, 1999-01-01 to 2023-12-31
2. Configure and Run LPU#
method = LPUMethod(n_members=200, station_id="02296750_peace_river")
result = method.run(obs, model, n_members=200, seed=42)
print(f"LPU diagnostics:")
print(f" fitted a : {result.config['fitted_a']:.6f}")
print(f" fitted b : {result.config['fitted_b']:.6f} CMS")
print(f" sigma2 : {result.config['sigma2']:.6f}")
print(f" kappa(JtJ) : {result.config['kappa']:.3e}")
print(f" SVD fallback : {result.config['svd_fallback']}")
print(f"Ensemble shape : {result.members.shape}")
LPU diagnostics:
fitted a : 0.893994
fitted b : -0.062006 CMS
sigma2 : 358.190624
kappa(JtJ) : 5.091e+03
SVD fallback : False
Ensemble shape : (200, 9131)
3. CI Band Visualization#
5th-95th percentile shaded band from the 200-member LPU ensemble.
p5 = np.percentile(result.members, 5, axis=0)
p95 = np.percentile(result.members, 95, axis=0)
fig, ax = plt.subplots(figsize=(14, 5))
ax.fill_between(result.time_index, p5, p95, alpha=0.3, color="blue", label="5th-95th CI")
ax.plot(result.time_index, result.model_output.values, color="green", linewidth=1.0, label="WAM")
ax.scatter(result.time_index, result.observations.values, s=6, color="red", alpha=0.5, label="Obs", zorder=5)
ax.set_xlabel("Date")
ax.set_ylabel("Discharge (CMS)")
ax.set_title("LPU Ensemble -- 5th to 95th Percentile CI Band")
ax.legend()
fig.tight_layout()
plt.show()
4. Export CSVs#
output_dir = Path("data/outputs")
ens_path = write_ensemble_csv(result, output_dir / "ensembles")
stats_path = write_stats_csv(result, output_dir / "stats")
print(f"Ensemble CSV: {ens_path}")
print(f"Stats CSV: {stats_path}")
Ensemble CSV: data\outputs\ensembles\lpu_02296750_peace_river_discharge_members.csv
Stats CSV: data\outputs\stats\lpu_02296750_peace_river_discharge_stats.csv
Summary#
Forward model: Q_sim = a * WAM + b (same as GLUE, different uncertainty quantification)
Fitting: scipy.optimize.leastsq with full_output=True
Parameter covariance: pcov = sigma2 * inv(JtJ)
Condition number check: kappa(JtJ) > 1e8 triggers SVD pseudo-inverse fallback
Ensemble: 200 members via MVN sampling of [a, b] parameter pairs
Output files:
data/outputs/ensembles/anddata/outputs/stats/