Linear Regression with PEST++ iES

Test notebook to compare the PEST++ iES outcomes for a linear regression problem where the confidence intervals can be computed exact.

Setup

Packages

import json
from dataclasses import dataclass, field
from pathlib import Path
from typing import Any, Literal

import matplotlib as mpl
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pyemu
import scipy as sp
import xarray as xr
from cmcrameri import cm as cmc
from numpy.typing import NDArray
from pastas import show_versions
from pyemu.utils.get_pestpp import run_main as get_pestpp

from pastas_plugins import show_plugin_versions
from pastas_plugins.pest import PestIesSolver

print("Pyemu version:", pyemu.__version__)
show_versions()
show_plugin_versions()

/home/docs/checkouts/readthedocs.org/user_builds/pastas-plugins/envs/latest/lib/python3.11/site-packages/tqdm/auto.py:21: TqdmWarning: IProgress not found. Please update jupyter and ipywidgets. See https://ipywidgets.readthedocs.io/en/stable/user_install.html

Pyemu version: 1.4.0
Pastas version: 1.13.2
Python version: 3.11.12
NumPy version: 2.4.3
Pandas version: 2.3.3
SciPy version: 1.17.1
Matplotlib version: 3.10.8
Numba version: 0.64.0

---------------------------------------------------------------------------
ImportError                               Traceback (most recent call last)
Cell In[1], line 23
     21 print("Pyemu version:", pyemu.__version__)
     22 show_versions()
---> 23 show_plugin_versions()

File ~/checkouts/readthedocs.org/user_builds/pastas-plugins/envs/latest/lib/python3.11/site-packages/pastas_plugins/__init__.py:27, in show_plugin_versions()
     25 for plugin in plugins:
     26     try:
---> 27         module = import_module(f"pastas_plugins.{plugin}.version")
     28         version = module.__version__
     29     except ModuleNotFoundError:

File ~/.asdf/installs/python/3.11.12/lib/python3.11/importlib/__init__.py:126, in import_module(name, package)
    124             break
    125         level += 1
--> 126 return _bootstrap._gcd_import(name[level:], package, level)

File <frozen importlib._bootstrap>:1204, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:1176, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:1126, in _find_and_load_unlocked(name, import_)

File <frozen importlib._bootstrap>:241, in _call_with_frames_removed(f, *args, **kwds)

File <frozen importlib._bootstrap>:1204, in _gcd_import(name, package, level)

File <frozen importlib._bootstrap>:1176, in _find_and_load(name, import_)

File <frozen importlib._bootstrap>:1147, in _find_and_load_unlocked(name, import_)

File <frozen importlib._bootstrap>:690, in _load_unlocked(spec)

File <frozen importlib._bootstrap_external>:940, in exec_module(self, module)

File <frozen importlib._bootstrap>:241, in _call_with_frames_removed(f, *args, **kwds)

File ~/checkouts/readthedocs.org/user_builds/pastas-plugins/envs/latest/lib/python3.11/site-packages/pastas_plugins/modflow/__init__.py:9
      1 # ruff: noqa: F401
      2 from pastas_plugins.modflow.modflow import (
      3     ModflowDrn,
      4     ModflowDrnSto,
   (...)      7     ModflowUzf,
      8 )
----> 9 from pastas_plugins.modflow.stressmodels import ModflowModel
     10 from pastas_plugins.modflow.version import __version__

File ~/checkouts/readthedocs.org/user_builds/pastas-plugins/envs/latest/lib/python3.11/site-packages/pastas_plugins/modflow/stressmodels.py:7
      5 from pastas.stressmodels import StressModelBase
      6 from pastas.timeseries import TimeSeries
----> 7 from pastas.typing import ArrayLike, TimestampType
      9 from .modflow import Modflow
     11 logger = getLogger(__name__)

ImportError: cannot import name 'TimestampType' from 'pastas.typing' (/home/docs/checkouts/readthedocs.org/user_builds/pastas-plugins/envs/latest/lib/python3.11/site-packages/pastas/typing/__init__.py)

Linear regression model

Setup of the model might be a bit weird but this way we can easily use the PestIesSolver in Pastas Plugins. The PestIesSolver allows for the pypestworker etc.

@dataclass
class LinRegModel:
    y: pd.Series = field(repr=False)

    def __post_init__(self):
        self.x: NDArray[float] = self.y.index.values
        p0_mean = np.mean(self.y.values)
        p1_ini = (y.iat[-1] - y.iat[0]) / (x[-1] - x[0])  # initial estimate slope
        self._parameters: pd.DataFrame = pd.DataFrame(
            [
                [
                    p0_mean,
                    np.min(self.y.values) - p0_mean / 4,
                    np.max(self.y.values) + p0_mean / 4,
                    np.nan,
                    True,
                ],
                [p1_ini, p1_ini - p1_ini * 2.0, p1_ini + p1_ini * 2.0, np.nan, True],
            ],
            index=pd.Index(
                [
                    "p_0",
                    "p_1",
                ],
                name="parnames",
            ),
            columns=["initial", "pmin", "pmax", "optimal", "vary"],
        )

    @property
    def p_0(self):
        return (
            self._parameters.at["p_0", "initial"]
            if np.isnan(self._parameters.at["p_0", "optimal"])
            else self._parameters.at["p_0", "optimal"]
        )

    @property
    def p_1(self):
        return (
            self._parameters.at["p_1", "initial"]
            if np.isnan(self._parameters.at["p_1", "optimal"])
            else self._parameters.at["p_1", "optimal"]
        )

    @property
    def parameters(self):
        return self._parameters

    def observations(self):
        return self.y

    def simulate(self):
        y = self.p_0 + self.p_1 * self.x
        return pd.Series(y, index=pd.Index(self.x, name="x"), name="Simulated")

    @staticmethod
    def synthetic(
        x: NDArray[float],
        p_0: float,
        p_1: float,
        noise_std: float = 0.0,
        noise_ccoeff: float = 0.0,
    ):
        """Generate synthetic linear data with optional noise.

        Parameters
        ----------
        x : NDArray[float]
            Input x values.
        p_0 : float
            Intercept of the linear model.
        p_1 : float
            Slope of the linear model.
        noise_std : float, optional
            Standard deviation of the Gaussian noise to be added, by default 0.0.
        noise_ccoeff : float, optional
            Correlation coefficient for autocorrelated noise (between -1 and 1), by default 0.0.

        Returns
        -------
        pd.Series
            Synthetic y values with optional noise.
        """
        y = p_0 + p_1 * x
        if noise_std > 0.0:
            drng = np.random.default_rng(pyemu.en.SEED)  # set seed
            noise = drng.normal(loc=0.0, scale=noise_std, size=len(y))
            if noise_ccoeff != 0.0:
                sige = np.sqrt(1 - noise_ccoeff**2) * noise_std
                e = drng.normal(loc=0.0, scale=sige, size=len(y))
                for j in range(1, len(y)):
                    noise[j] = noise_ccoeff * noise[j - 1] + e[j]

            y += noise
        return pd.Series(y, index=pd.Index(x, name="x"))

    def set_parameter(
        self,
        name: str,
        initial: float | None = None,
        pmin: float | None = None,
        pmax: float | None = None,
        optimal: float | None = None,
        vary: bool | None = None,
    ) -> None:
        if initial is not None:
            self._parameters.loc[name, "initial"] = initial
        if pmin is not None:
            self._parameters.loc[name, "pmin"] = pmin
        if pmax is not None:
            self._parameters.loc[name, "pmax"] = pmax
        if optimal is not None:
            self._parameters.loc[name, "optimal"] = optimal
            # setattr(self, name, optimal)
        if vary is not None:
            self._parameters.loc[name, "vary"] = vary

    def add_solver(self, solver: Any):
        self.solver = solver
        solver.set_model(self)

    def solve(self) -> None:
        self.solver.solve()

    def to_file(self, path: Path) -> None:
        # write dummy model to emulate pastas
        with open(path, "w") as f:
            f.write("Dummy model file\n")
            f.write(f"p_0: {self.p_0}\n")
            f.write(f"p_1: {self.p_1}\n")

    def rss(self, sim: pd.Series | None = None) -> float:
        if sim is None:
            sim = self.simulate()
        residuals = self.y - sim
        return float(np.sum(residuals**2))

Linear regression solver

Linear regression solver from SciPy (scipy.stats.linregress) which can compute the true confidence intervals and everything.

@dataclass
class LinRegSolver:
    """https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.linregress.html"""

    model: None | LinRegModel = None

    def set_model(self, model: LinRegModel) -> None:
        self.model = model

    @staticmethod
    def linregress(x: NDArray[float], y: NDArray[float]) -> sp.optimize.OptimizeResult:
        res = sp.stats.linregress(x=x, y=y)
        return res

    def solve(self) -> None:
        if self.model is None:
            raise ValueError("No model assigned to solver.")
        res = self.linregress(x=self.model.x, y=self.model.y)
        self.res = res
        self.model.set_parameter(name="p_1", optimal=res.slope)
        self.model.set_parameter(name="p_0", optimal=res.intercept)

    def critical_t(self, alpha: float = 0.05) -> float:
        """Calculate the critical t-value for the slope and intercept.

        Parameters
        ----------
        alpha : float, optional
            Significance level, by default 0.05

        Returns
        -------
        float
            Critical t-value
        """
        n = len(self.model.x)
        dof = n - 2  # degrees of freedom
        # alpha = 0.05 # probability of confidence interval
        tinv = abs(
            sp.stats.distributions.t.ppf(alpha / 2.0, dof)
        )  # Two-sided inverse Students t-distribution
        return tinv

    def ci(self, alpha: float = 0.05) -> tuple[float, float]:
        """Calculate the confidence interval for the slope and intercept.

        Parameters
        ----------
        alpha : float, optional
            Significance level, by default 0.05

        Returns
        -------
        tuple[float, float]
            Confidence interval for intercept and slope (p_0_ci, p_1_ci)
        """
        tinv = self.critical_t(alpha=alpha)
        p_0_ci = tinv * self.res.intercept_stderr
        p_1_ci = tinv * self.res.stderr
        return p_0_ci, p_1_ci

    def ci_sim(self, alpha: float = 0.05) -> tuple[NDArray[float], NDArray[float]]:
        """Calculate the confidence interval for the simulated values.
        https://rpubs.com/aaronsc32/regression-confidence-prediction-intervals

        Parameters
        ----------
        alpha : float, optional
            Significance level, by default 0.05

        Returns
        -------
        tuple[NDArray[float], NDArray[float]]
            Lower and upper confidence interval for simulated values
        """
        x = self.model.x
        n = len(x)
        se = np.sqrt(self.model.rss() / (n - 2)) * np.sqrt(
            1 / n + (x - np.mean(x)) ** 2 / np.sum((x - np.mean(x)) ** 2)
        )

        y_fit = self.model.simulate()
        tinv = self.critical_t(alpha=alpha)
        lb = y_fit - tinv * se
        ub = y_fit + tinv * se

        # Compute lower and upper bounds for simulated values using the confidence intervals
        # p0, p1 = self.res.intercept, self.res.slope
        # p_0_ci, p_1_ci = self.ci(alpha=alpha)
        # b0 = self.model.synthetic(self.model.x, p0 - p_0_ci, p1 - p_1_ci).rename("--")
        # b1 = self.model.synthetic(self.model.x, p0 + p_0_ci, p1 + p_1_ci).rename("++")
        # b2 = self.model.synthetic(self.model.x, p0 - p_0_ci, p1 + p_1_ci).rename("-+")
        # b3 = self.model.synthetic(self.model.x, p0 + p_0_ci, p1 - p_1_ci).rename("+-")
        # lb = pd.concat([b0, b1, b2, b3], axis=1).min(axis=1)
        # ub = pd.concat([b0, b1, b2, b3], axis=1).max(axis=1)

        return lb, ub

PEST++ iES PostProcessor

Load all iES results in one object (xarray.Dataset).

class PestIesPostProcessor:
    """PEST++ iES postprocessing class"""

    def __init__(self, path: str | Path):
        """Initialize the PestIesPostProcessor with a path."""
        self.path = Path(path).resolve()

    def __repr__(self):
        return f"PestIesPostProcessor(path={self.path})"

    @property
    def files(self) -> list[Path]:
        """Return all files in the path."""
        return [x for x in self.path.glob("*") if x.is_file()]

    @property
    def noptmax(self) -> int:
        """Return the maximum number of optimization iterations."""
        nopt_maxes = [x.name.split(".")[1] for x in list(self.path.glob("pest.*.*"))]
        noptmax = max([int(x) for x in nopt_maxes if x.isdigit()])
        return noptmax

    @staticmethod
    def df_to_da(df: pd.DataFrame) -> xr.DataArray:
        """Convert a DataFrame to an xarray DataArray."""
        da = xr.DataArray(
            df.values,
            coords=[df.index, df.columns],
            dims=[
                "index" if df.index.name is None else df.index.name,
                "columns" if df.columns.name is None else df.columns.name,
            ],
        )
        return da

    def _load_ensembles(self, suffix: Literal[".par.csv", ".obs.csv"]) -> xr.Dataset:
        """Load ensembles from the processed files."""
        suffix_mapping = {
            ".obs.csv": "oname",
            ".par.csv": "pname",
            ".pdc.csv": "pdcname",
        }

        das = []
        for file in self.files:
            if file.name.endswith(suffix):
                if "rejected" in file.name:
                    continue
                df = pd.read_csv(file, index_col=0)
                df.columns.name = suffix_mapping[suffix]
                iteration = int(file.name.split(".")[1])
                da = self.df_to_da(df).rename(suffix.split(".")[1])
                da = da.expand_dims(iteration=[iteration])
                das.append(da)
        return xr.concat(das, dim="iteration").sortby("iteration")

    def _load_file(self, file: Path, column_name: str) -> pd.DataFrame:
        if file not in self.files:
            raise FileNotFoundError(f"File not found: {file} in {self.path}")
        df = pd.read_csv(file, index_col=0)
        df.columns.name = column_name
        return df

    def load_parameters(self) -> xr.Dataset:
        """Load parameter ensembles from the processed files."""
        return self._load_ensembles(".par.csv")

    def load_simulations(self) -> xr.Dataset:
        """Load simulation ensembles from the processed files."""
        return self._load_ensembles(".obs.csv")

    def load_observations(self) -> xr.DataArray:
        df = self._load_file(self.path / "pest.obs+noise.csv", "oname")
        da = self.df_to_da(df).rename("observations")
        return da

    def load_weights(self) -> xr.DataArray:
        """Load weights from the processed files."""
        df = self._load_file(self.path / "pest.weights.csv", "oname")
        df.index = df.index.astype(str)
        df = df.rename(index={df.index[-1]: "base"})
        da = self.df_to_da(df).rename("weights")
        return da

    def load_phi(self) -> xr.DataArray:
        """Load phi ensembles from the processed files."""
        df = self._load_file(self.path / "pest.phi.actual.csv", "real_name")
        col_i = df.columns.to_list().index("max")  # start reading after the max column
        df = df.iloc[:, slice(col_i + 1, None)]
        da = self.df_to_da(df).rename("phi")
        return da

    def load_dataset(self, rename: bool = False) -> xr.Dataset:
        """Load the complete dataset from the processed files."""
        ds = xr.merge(
            [
                self.load_parameters(),
                self.load_simulations(),
                self.load_phi(),
                self.load_observations(),
                self.load_weights(),
            ],
            compat="override",
        )
        if rename:
            param_index = json.load((self.path / "parameter_index.json").open("r"))
            obs_index = json.load((self.path / "observation_index.json").open("r"))
            ds = ds.assign_coords(
                pname=list(param_index.values()),
                oname=[pd.Timestamp(x) for x in obs_index.values()],
            )
        return ds

Case 1: no noise on observations

Synthetic series with no noise

# create synthetic data
x = np.linspace(0.0, 1.0, 101)  # x coordinates
a = 1.0  # slope
b = 5.0  # y-intercept
noise_std = 0.0  # noise level

y = LinRegModel.synthetic(p_0=b, p_1=a, x=x, noise_std=noise_std)

Solve model

SciPy

# create a linear regression model and solve with SciPy
ml_sp = LinRegModel(y)
ml_sp.add_solver(LinRegSolver())
ml_sp.solve()
ml_sp.parameters

	initial	pmin	pmax	optimal	vary
parnames
p_0	5.5	3.625	7.375	5.0	True
p_1	1.0	-1.000	3.000	1.0	True

Scipy finds the true parameters.

Results

Let’s look at the performance of both models.

Simulated values

fig, ax = plt.subplots(figsize=(6.75, 3.0))
ax.plot(
    y.index,
    ds["obs"].isel(iteration=-1).values.T,
    color="C0",
    alpha=0.1,
    label="iES",
)
ax.plot(x, y, label="Observations", color="k", marker=".", linestyle="none")
ax.plot(
    x, ml_sp.simulate(), color="C1", linestyle="--", label="SciPy Linear Regression"
)
lower_bound, upper_bound = ml_sp.solver.ci_sim()
ax.fill_between(
    x,
    lower_bound,
    upper_bound,
    color="C1",
    alpha=0.2,
    label=f"{0.95:0.0%} Confidence Interval",
)
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles[-4:], labels[-4:])
ax.grid()

../_images/e2aa2c4949c41b84b4d056f24d4d5e173f586febbafdc65985859a4cf085472c.png

SciPy is perfect on the black observations, confidence interval is the same . PEST++ iES is much wider but around the black observations.

Parameter estimates

mosaic = [["p_0"], ["p_1"]]
f, axd = plt.subplot_mosaic(
    mosaic,
    figsize=(7.0, 6.0),
    sharex=True,
    layout="tight",
)
scipy_ci = ml_sp.solver.ci()
for i, pname in enumerate(axd):
    par = ds["par"].sel(pname=pname)
    parq = par.quantile(
        [0.025, 0.975], dim="real_name", skipna=True
    )  # .values.reshape(2, -1)
    axd[pname].fill_between(
        parq["iteration"],
        parq.values[0],
        parq.values[1],
        color="C0",
        label="iES 95% CI",
        linewidth=2.0,
    )

    scipy_opt = ml_sp.parameters.at[pname, "optimal"]
    axd[pname].axhline(
        scipy_opt,
        color="C1",
        alpha=1.0,
        linewidth=2.0,
        label="SciPy 95% CI\n(no iterations)",
    )
    axd[pname].fill_between(
        x=[0.0, ml_ies_pp.noptmax],
        y1=[scipy_opt - scipy_ci[i], scipy_opt - scipy_ci[i]],
        y2=[scipy_opt + scipy_ci[i], scipy_opt + scipy_ci[i]],
        color="C1",
        alpha=0.2,
    )
    axd[pname].set(
        xlabel="Iteration" if "p_1" in pname else None,
        ylabel=pname,
        xlim=(0, ml_ies_pp.noptmax),
        ylim=parq.sel(iteration=1),
    )
    axd[pname].xaxis.set_major_locator(mpl.ticker.MultipleLocator(5.0))
    axd[pname].xaxis.set_minor_locator(mpl.ticker.MultipleLocator(1.0))
axd["p_0"].legend()

<matplotlib.legend.Legend at 0x753addbb2710>

../_images/3928632bf531f1e69c9105c673806eaafb740d4cf49decaac87cdb690e079ab3.png

iES stops after 21 iterations

Response surface

def get_response_surface(ml: LinRegModel, N: int) -> NDArray[float]:
    """Calculate the response surface for parameters p_0 and p_1."""
    p0 = np.linspace(
        ml.parameters.at["p_0", "pmin"], ml.parameters.at["p_0", "pmax"], N
    )
    p1 = np.linspace(
        ml.parameters.at["p_1", "pmin"], ml.parameters.at["p_1", "pmax"], N
    )
    P0, P1 = np.meshgrid(p0, p1)
    RSS = np.zeros(P0.shape)
    for i in range(N):
        for j in range(N):
            sim = ml.synthetic(ml.x, P0[i, j], P1[i, j])
            RSS[i, j] = ml.rss(sim)
    return P0, P1, RSS


def plot_response_surface(
    P0: NDArray[float],
    P1: NDArray[float],
    RSS: NDArray[float],
    ax: plt.Axes,
    vmin: float | None = None,
    vmax: float | None = None,
) -> plt.Axes:
    """Plot the response surface for parameters p_0 and p_1."""
    if vmin is None:
        vmin = np.min(RSS)
    if vmax is None:
        vmax = np.max(RSS)

    exp_min = int(np.floor(np.log10(vmin)))  # round down to nearest integer
    exp_max = int(np.ceil(np.log10(vmax)))  # round up to nearest integer

    levels = np.logspace(exp_min, exp_max, num=(exp_max - exp_min + 1))
    norm = mpl.colors.LogNorm(vmin=vmin, vmax=vmax)
    cs = ax.pcolormesh(P0, P1, RSS, cmap=cmc.batlow, norm=norm, shading="nearest")
    ct = ax.contour(P0, P1, RSS, levels=levels, colors="w")
    ax.clabel(ct, levels, fmt="%.1e")

    ax.get_figure().colorbar(cs, ax=ax, label="Sum of Squared Residuals")
    ax.set_xlabel("p_0")
    ax.set_ylabel("p_1")
    ax.set_title("Response Surface", loc="left")
    return ax

P0, P1, RSS = get_response_surface(ml_sp, 500)

f, ax = plt.subplots(figsize=(6.75, 5.0))
ax.set_xlim(4.8, 5.2)
ax.set_ylim(0.6, 1.4)
plot_response_surface(P0, P1, RSS, vmin=1e-3, vmax=1e1, ax=ax)
ax.axhline(a, color="k", linewidth=0.5, linestyle="--", label="True")
ax.axvline(b, color="k", linewidth=0.5, linestyle="--")

traj_kwargs = {
    "color": "w",
    "marker": ".",
    "linewidth": 0.5,
    "markersize": 2.0,
    "alpha": 0.2,
}
# for real in ds["par"].real_name.values:
#     ax.plot(
#         ds["par"].sel(pname="p_0", real_name=real),
#         ds["par"].sel(pname="p_1", real_name=real),
#         **traj_kwargs,
#     )
ax.plot(
    ds["par"].isel(iteration=-1).sel(pname="p_0"),
    ds["par"].isel(iteration=-1).sel(pname="p_1"),
    **traj_kwargs,
    linestyle="none",
)

sp_ci = ml_sp.solver.ci()
ax.errorbar(
    x=ml_sp.p_0,
    y=ml_sp.p_1,
    xerr=sp_ci[0],
    yerr=sp_ci[1],
    linewidth=0.8,
    marker="o",
    markersize=1.0,
    color="C1",
    label="SciPy 95% CI",
)

ies_ci = ds["par"].isel(iteration=-1).quantile([0.025, 0.975], "real_name").values
ies_base = ds["par"].isel(iteration=-1).sel(real_name="base").values
ax.errorbar(
    x=ies_base[0],
    y=ies_base[1],
    xerr=np.abs(ies_ci[:, [0]] - ies_base[0]),
    yerr=np.abs(ies_ci[:, [1]] - ies_base[1]),
    color="C0",
    linewidth=0.8,
    label="iES 95% CI",
)
ax.legend(loc="lower right", bbox_to_anchor=(1.0, 1.0), frameon=False)

<matplotlib.legend.Legend at 0x753add247890>

../_images/c5a5654ccd2fbb2730f48d239cc856c6b43cdfabde44ce3e6ed615e1039a0b5c.png

Loss function

f, ax = plt.subplots(figsize=(6.75, 3.0))
ax.plot(ds["phi"].mean("real_name") / weight**2, label="iES mean")
ax.fill_between(
    ds["phi"].iteration,
    ds["phi"].min("real_name") / weight**2,
    ds["phi"].max("real_name") / weight**2,
    color="C0",
    alpha=0.2,
    label="iES min/max bounds",
)
print(f"SciPy Phi: {ml_sp.rss():.2e}")
# ax.axhline(ml_sp.rss(), color="C1", linestyle="--", label="SciPy RSS")
ax.plot([], [], color="C1", linestyle="--", label=f"SciPy: {ml_sp.rss():0.2e}")
ax.semilogy()
ax.set(
    xlim=(0, ml_ies_pp.noptmax),
    # ylim=(1e-4, 1e2),
    ylabel="Phi | RSS",
    xlabel="iteration",
)
ax.xaxis.set_major_locator(mpl.ticker.MultipleLocator(5.0))
ax.xaxis.set_minor_locator(mpl.ticker.MultipleLocator(1.0))
ax.legend()

SciPy Phi: 2.00e-28

<matplotlib.legend.Legend at 0x753add11c550>

../_images/34683bcc141a67eccfd45251570a18dd44d67478564bac24fb59e3cfa387c7be.png

The fit for both iES is still a lot larger than SciPy but the outcomes are pretty much the same. The CI is basically zero for both.

Case 2: noise on observations

noise_std = 0.05  # noise level
y_n = LinRegModel.synthetic(p_0=b, p_1=a, x=x, noise_std=noise_std)

Solve model

SciPy

ml_sp_n = LinRegModel(y_n)
ml_sp_n.add_solver(LinRegSolver())
ml_sp_n.solve()
ml_sp_n.parameters

	initial	pmin	pmax	optimal	vary
parnames
p_0	5.501794	3.63682	7.386489	5.015022	True
p_1	1.000000	-1.00000	3.000000	0.973545	True

Results

Simulated values

fig, ax = plt.subplots(figsize=(6.75, 4.0))
ax.plot(
    y_n.index,
    ds_n["obs"].isel(iteration=-1).values.T,
    color="C0",
    alpha=0.1,
    label="iES",
)
ax.plot(x, y_n, label="Observations", color="k", marker=".", linestyle="none")
ax.plot(
    x, ml_sp_n.simulate(), color="C1", linestyle="--", label="SciPy Linear Regression"
)
lower_bound, upper_bound = ml_sp_n.solver.ci_sim()
ax.fill_between(
    x,
    lower_bound,
    upper_bound,
    color="C1",
    alpha=0.3,
    label=f"{0.95:0.0%} Confidence Interval",
    zorder=10,
)
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles[-4:], labels[-4:])
ax.grid(True)

../_images/1465536cbf7432398a99ae18a1a2e283fdf140e568553b183474b887afc4faf9.png

Parameter estimates

mosaic = [["p_0"], ["p_1"]]
f, axd = plt.subplot_mosaic(
    mosaic,
    figsize=(7.0, 6.0),
    sharex=True,
    layout="tight",
)
scipy_ci = ml_sp_n.solver.ci()
for i, pname in enumerate(axd):
    par = ds_n["par"].sel(pname=pname)
    parq = par.quantile(
        [0.025, 0.975], dim="real_name", skipna=True
    )  # .values.reshape(2, -1)
    axd[pname].fill_between(
        parq["iteration"],
        parq.values[0],
        parq.values[1],
        color="C0",
        label="iES 95% CI",
        linewidth=2.0,
    )

    scipy_opt = ml_sp_n.parameters.at[pname, "optimal"]
    axd[pname].axhline(
        scipy_opt,
        color="C1",
        alpha=1.0,
        linewidth=2.0,
        label="SciPy 95% CI\n(no iterations)",
    )
    axd[pname].fill_between(
        x=[0.0, ml_ies_pp_n.noptmax],
        y1=[scipy_opt - scipy_ci[i], scipy_opt - scipy_ci[i]],
        y2=[scipy_opt + scipy_ci[i], scipy_opt + scipy_ci[i]],
        color="C1",
        alpha=0.5,
    )
    axd[pname].xaxis.set_major_locator(mpl.ticker.MultipleLocator(10.0))
    axd[pname].xaxis.set_minor_locator(mpl.ticker.MultipleLocator(1.0))
    axd[pname].set_ylim(parq.sel(iteration=1))
    axd[pname].set_ylabel(pname)
    axd[pname].set_xlabel("Iteration") if "p_1" in pname else None
    axd[pname].set_xlim(0, ml_ies_pp_n.noptmax)
axd["p_0"].legend()

<matplotlib.legend.Legend at 0x753adcebed50>

../_images/12a07c4fe84ea12b5b1b6e79654ec5f0f311ba31652e5be604660afc75363830.png

Response surface

P0_n, P1_n, RSS_n = get_response_surface(ml_sp_n, 500)

f, ax = plt.subplots(figsize=(6.75, 5.0))
ax.set(
    xlim=(4.8, 5.2),
    ylim=(0.6, 1.4),
)
plot_response_surface(P0_n, P1_n, RSS_n, vmin=np.min(RSS_n), vmax=1e1, ax=ax)
ax.axhline(a, color="k", linewidth=0.5, linestyle="--", label="True")
ax.axvline(b, color="k", linewidth=0.5, linestyle="--")

traj_kwargs = {
    "color": "w",
    "marker": ".",
    "linewidth": 0.5,
    "markersize": 2.0,
    "alpha": 0.2,
}
# for real in ds_n["par"].real_name.values:
#     ax.plot(
#         ds_n["par"].sel(pname="p_0", real_name=real),
#         ds_n["par"].sel(pname="p_1", real_name=real),
#         **traj_kwargs,
#     )
ax.plot(
    ds_n["par"].isel(iteration=-1).sel(pname="p_0"),
    ds_n["par"].isel(iteration=-1).sel(pname="p_1"),
    **traj_kwargs,
    linestyle="none",
)

sp_ci = ml_sp_n.solver.ci()
ax.errorbar(
    x=ml_sp_n.p_0,
    y=ml_sp_n.p_1,
    xerr=sp_ci[0],
    yerr=sp_ci[1],
    linewidth=0.8,
    color="C1",
    label="SciPy 95% CI",
)

ies_ci = ds_n["par"].isel(iteration=-1).quantile([0.025, 0.975], "real_name").values
ies_base = ds_n["par"].isel(iteration=-1).sel(real_name="base").values
ax.errorbar(
    x=ies_base[0],
    y=ies_base[1],
    xerr=np.abs(ies_ci[:, [0]] - ies_base[0]),
    yerr=np.abs(ies_ci[:, [1]] - ies_base[1]),
    color="C0",
    linewidth=0.8,
    label="iES 95% CI",
)
ax.legend(loc="lower right", bbox_to_anchor=(1.0, 1.0), frameon=False)

<matplotlib.legend.Legend at 0x753adcc11a90>

../_images/6cabd0c9ee7b3504ff4c11dcdb7dd4641acb74750b7c83fcaab405b6d2879a3b.png

Outcomes are pretty similar for iES and SciPy again. The confidence intervals look pretty much the same.

Loss function

Loss function in PEST++ iES

f, ax = plt.subplots(figsize=(6.75, 3.5))
ax.plot(ds_n["phi"].mean("real_name") / weight**2, label="iES mean")
ax.fill_between(
    ds_n["phi"].iteration,
    ds_n["phi"].min("real_name") / weight**2,
    ds_n["phi"].max("real_name") / weight**2,
    color="C0",
    alpha=0.2,
    label="iES min/max bounds",
)
ax.axhline(ml_sp_n.rss(), color="C1", linestyle="--", label="SciPy RSS")
ax.set(
    xlim=(0, ml_ies_pp_n.noptmax),
    # ylim=(2e-1, 1e0),
    ylabel="Phi | RSS",
    xlabel="iteration",
)
ax.semilogy()
ax.xaxis.set_major_locator(mpl.ticker.MultipleLocator(10.0))
ax.xaxis.set_minor_locator(mpl.ticker.MultipleLocator(1.0))

ax.legend()

<matplotlib.legend.Legend at 0x753adcba39d0>

../_images/f2f64a6095e5bd32dc28bb4f2a2da4ebc4794761ea9c3853fed8a78b8d67649f.png

Linear Regression with PEST++ iES

Setup

Packages

Linear regression model

Linear regression solver

PEST++ iES PostProcessor

Case 1: no noise on observations

Solve model

SciPy

PEST++ iES

Results

Simulated values

Parameter estimates

Response surface

Loss function

Case 2: noise on observations

Solve model

SciPy

PEST++ iES

Results

Simulated values

Parameter estimates

Response surface

Loss function

Case 3: noise on observations but not on iES realizations

Solve model

PEST++ iES

Results

Simulated values

Parameter estimates

Response surface

Loss function