A small estimation-to-design example - Optimization for Decision Science

One biochemical oxygen demand (BOD) model runs through parameter estimation, uncertainty analysis, and experimental design. The complete Pyomo model, calculations, and diagnostics remain together so that a later step uses the results of the earlier fit.

Part I: parameter estimation and local uncertainty. Specify the measurement model, derive the likelihood, fit the data, differentiate the predictions, and construct the local confidence ellipse.

Part II: nonlinear uncertainty and experimental design. Compare likelihood-ratio regions, bootstrap fits, and profiles; choose a new measurement time; then measure, append the data, and refit. The two larger reaction-kinetics tutorials apply the same workflow to dynamic models.

The historical data, the assumed measurement-error model, and the digital experiment serve different roles. The six historical observations determine the initial fit. The known standard deviation of 1 mg/L is an illustrative assumption. The digital experiment supplies new simulated measurements when testing a design.

import sys

if "google.colab" in sys.modules:
    !wget -q "https://raw.githubusercontent.com/ndcbe/optimization/main/notebooks/helper.py"
    import helper
    helper.easy_install()
else:
    sys.path.insert(0, "../")
    import helper

helper.set_plotting_style()

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pyomo.environ as pyo
from pyomo.contrib.doe import DesignOfExperiments, ObjectiveLib
from pyomo.contrib.parmest.experiment import Experiment
from pyomo.contrib.parmest.parmest import Estimator

Part I: parameter estimation and local uncertainty¶

Model and observations¶

The response model is

y(t;\boldsymbol{\theta})=\theta_1\left(1-\exp(-\theta_2t)\right),

(1)

where:

$\theta_1$ is the long-time asymptote;
$\theta_2$ is the approach rate; and
$t$ is the measurement time.

We assume independent measurement noise with known standard deviation $\sigma=1$ . In Pyomo’s experiment suffix, measurement_error means standard deviation, not variance or weight.

The legacy code column is named time_days; its values in this BOD example are days.

Time is measured in days, response and asymptote in mg/L, and rate in day $^{-1}$ . The column and Pyomo component are named time_days consistently. Laky et al. (2026), Section 4.1, Equation (32), gives the same model and units; its two-observation example is different from the six-observation fit here.

Source. This is Bates and Watts’ own BOD example. The six observations below are their Table A1.4 (Nonlinear Regression Analysis and Its Applications, Wiley, 1988, p. 270), and the response model above is their Equation (2.8), p. 41. The data are not originally theirs: they reproduce them from Marske (1967), Biochemical Oxygen Demand Data Interpretation Using Sum of Squares Surface, M.S. thesis, University of Wisconsin--Madison (their reference list, p. 339). Each entry is the average of two analyses on one inoculated bottle of stream water. Fitting these six points reproduces their published estimate $\hat{\boldsymbol{\theta}} = (19.143,\ 0.5311)$ with $s^2 = 6.498$ on four degrees of freedom (Example BOD 4, p. 51) -- which is what the assertion in the estimation cell checks. The same data set reaches this course through parmest’s Rooney--Biegler tutorial problem.

The response printed in Appendix A1.4, p. 270 omits the minus sign in the exponent; Equation (2.8), p. 41 has the decreasing exponential used here. A positive exponent with positive rate would give a negative, decreasing demand and is not the fitted BOD model.

data = pd.DataFrame(
    {
        "time_days": [1, 2, 3, 4, 5, 7],
        "y": [8.3, 10.3, 19.0, 16.0, 15.6, 19.8],
    }
)
data

A digital experiment¶

At the end of this notebook we will design a new measurement and then take it. The BOD data set is historical, so there is no seventh flask to go back to. We therefore build a digital experiment: a known ground truth $\boldsymbol{\theta}^{\text{true}}$ plus Gaussian noise, which can be queried at any time we like.

The ground truth is chosen close to the historical fit to illustrate the workflow. This choice does not validate the Gaussian noise model against the historical data.

Two data sets are now in play, and it is worth being blunt about which is which:

data -- the real BOD observations. This stays the working data set for every estimate below, so the published numbers do not move.
simulated -- six draws from the digital experiment at the same times. It appears once to illustrate the generator and is never used for estimation.

The digital experiment earns its keep in the final section, where it supplies the one new measurement the design step asks for.

# Ground truth for the digital experiment. Deliberately close to the historical
# fit; this illustrates the workflow without validating the noise assumption.
THETA_TRUE = {"asymptote": 19.0, "rate_constant": 0.53}
MEASUREMENT_STD = 1.0


def bod_response(t, theta):
    """Noise-free response y(t) = theta_1 * (1 - exp(-theta_2 * t))."""
    t = np.asarray(t, dtype=float)
    return theta["asymptote"] * (1.0 - np.exp(-theta["rate_constant"] * t))


def run_digital_experiment(t, rng):
    """Take one measurement at time t: the true response plus one noise draw."""
    return float(bod_response(t, THETA_TRUE) + rng.normal(0.0, MEASUREMENT_STD))


# A fixed seed makes every number in this notebook reproducible, including the
# new measurement taken after the design step.
calibration_rng = np.random.default_rng(60499)

simulated = pd.DataFrame(
    {
        "time_days": data["time_days"],
        "y": [run_digital_experiment(t, calibration_rng) for t in data["time_days"]],
    }
)

# Side by side: the historical measurements and six draws from the simulator.
# Compare them without treating visual similarity as model validation.
data.rename(columns={"y": "y (real)"}).assign(**{"y (digital)": simulated["y"].round(2)})

From a measurement model to weighted least squares¶

For experiment $i$ , hold the design $\boldsymbol{\varphi}_i$ fixed and predict $m_i$ outputs:

\hat{\mathbf y}_i=\mathbf f(\boldsymbol{\varphi}_i,\boldsymbol\theta),\qquad \mathbf Y_i=\hat{\mathbf y}_i(\boldsymbol\theta)+\boldsymbol\varepsilon_i,\qquad \boldsymbol\varepsilon_i\sim\mathcal N(\mathbf0,\boldsymbol\Sigma_i).

(2)

The parameter $\boldsymbol\theta\in\Theta\subseteq\mathbb R^p$ is fixed but unknown. The estimator $\hat{\boldsymbol\theta}(\mathbf Y)$ is a rule applied to random data; its value for the observations $\mathbf y$ is an estimate. Assume independent experiments and known positive-definite $\boldsymbol\Sigma_i$ that do not depend on the parameters. Outputs within an experiment may be correlated.

Define $\mathbf r_i(\boldsymbol\theta)=\mathbf y_i-\hat{\mathbf y}_i(\boldsymbol\theta)$ . The likelihood is a function of the parameters with the observed data held fixed:

L(\boldsymbol\theta;\mathbf y) =\prod_{i=1}^N \frac{\exp[-\tfrac12\mathbf r_i^T\boldsymbol\Sigma_i^{-1}\mathbf r_i]} {(2\pi)^{m_i/2}(\det\boldsymbol\Sigma_i)^{1/2}}.

(3)

Taking the negative logarithm turns the product into a sum:

-\log L =\frac12\sum_i m_i\log(2\pi) +\frac12\sum_i\log\det\boldsymbol\Sigma_i +\frac12\sum_i\mathbf r_i^T\boldsymbol\Sigma_i^{-1}\mathbf r_i.

(4)

The first two terms are constant in $\boldsymbol\theta$ . Thus the maximum likelihood estimate minimizes

\ell(\boldsymbol\theta)=\frac12\sum_i \mathbf r_i^T\boldsymbol\Sigma_i^{-1}\mathbf r_i.

(5)

Here $\ell$ always means the summed negative log likelihood with parameter-independent constants removed. For diagonal covariance, this is $\tfrac12\sum_{i,j}(r_{ij}/\sigma_{ij})^2$ . Halving a standard deviation multiplies its weight by four. For the six BOD observations, $m_i=1$ and $\sigma=1$ mg/L, so

\ell(\boldsymbol\theta)=\frac12\sum_{i=1}^6 [y_i-\theta_1(1-e^{-\theta_2t_i})]^2/\sigma^2.

(6)

See Seber & Wild (1989), §2.1.4, p. 27 and §2.2.1, Eqs. (2.45)--(2.46), pp. 32--33; Nocedal & Wright (2006), §10.1, pp. 249--250. The number of Gaussian scalar components is $\sum_i m_i$ ; the number of experiments is $N$ . These counts coincide in this scalar-output example.

The estimation/design workflow: update the model or carry remaining uncertainty into a decision.

Generate measurements at fixed designs from a shared response model and independent experiment errors.

Known covariance standardizes the residual before its squared length is computed.

Change the likelihood when the error model changes¶

Stack all residuals in $\mathbf r$ . If errors across experiments are correlated, use their joint Gaussian covariance $\boldsymbol\Omega\succ0$ :

-\log L=\frac{n_y}{2}\log(2\pi)+\frac12\log\det\boldsymbol\Omega +\frac12\mathbf r^T\boldsymbol\Omega^{-1}\mathbf r, \qquad n_y=\sum_i m_i.

(7)

Off-diagonal blocks describe dependence between experiments; multiplying their marginal densities would discard this dependence. A scalar measurement_error suffix cannot encode a full off-diagonal covariance. Implement the joint likelihood explicitly for such a model.

If $\boldsymbol\Omega$ depends on an unknown parameter, retain both its log determinant and its inverse in the objective. Minimizing only the weighted residual term can favor artificially large variances. For independent scalar Laplace errors of known scale $b>0$ , the density is $\exp(-|r_i|/b)/(2b)$ , and the parameter-dependent objective is instead $\sum_i|r_i|/b$ .

If the response model is implicit, share one parameter vector across all experiments and fit states and predictions together:

\begin{aligned} \min_{\boldsymbol\theta,\{\mathbf x_i,\hat{\mathbf y}_i\}} &\quad\frac12\sum_i(\mathbf y_i-\hat{\mathbf y}_i)^T \boldsymbol\Sigma_i^{-1}(\mathbf y_i-\hat{\mathbf y}_i)\\ \text{subject to}&\quad \mathbf F_i(\mathbf x_i,\boldsymbol\varphi_i,\boldsymbol\theta,\hat{\mathbf y}_i)=\mathbf0, \quad i=1,\ldots,N,\qquad\boldsymbol\theta\in\Theta. \end{aligned}

(8)

A discretized DAE changes the model constraints, while the likelihood still follows from the measurement assumptions. The dynamic ParmEst tutorial uses this pattern. For other error laws, start with the density before choosing a residual norm; see Nocedal & Wright (2006), §10.1 and Seber & Wild (1989), Chapter 2.

Label one experiment¶

Pyomo’s Experiment abstraction uses four local suffixes:

Suffix	Role in this example
`experiment_outputs`	measured response and observed value
`unknown_parameters`	parameters and nominal values
`measurement_error`	response standard deviation
`experiment_inputs`	candidate measurement time

parmest uses the first three. Pyomo.DoE also needs experiment_inputs.

class BODExperiment(Experiment):
    """One biochemical oxygen demand measurement and its model labels.

    Parameters
    ----------
    data : mapping
        Measurement with keys ``time_days`` (day) and ``y`` (mg/L).
    measurement_standard_deviation : float, optional
        Known measurement standard deviation in mg/L; default 1 is illustrative.
    theta : mapping, optional
        Initial or nominal asymptote (mg/L) and rate_constant (1/day).
    design_bounds : tuple, optional
        Lower and upper allowed measurement times, in days.
    """

    def __init__(self, data, measurement_standard_deviation=1.0,
                 theta=None, design_bounds=(0, 10)):
        # Store the measurement, assumed noise, and initial parameter guesses.
        self.data = dict(data)
        self.sigma = measurement_standard_deviation
        self.theta = theta or {"asymptote": 15.0, "rate_constant": 0.5}
        self.design_bounds = design_bounds
        self.model = None

    def get_labeled_model(self):
        """Build and return the cached Pyomo model with four local suffixes."""
        if self.model is not None:
            return self.model
        # Construct the equations, then attach the reusable tool labels.
        m = build_bod_model(self.data, self.theta, self.design_bounds)
        self.model = label_bod_model(m, self.data, self.sigma)
        return self.model

def build_bod_model(data, theta, design_bounds):
    """Construct one response equation with physically consistent units."""
    m = pyo.ConcreteModel()
    # Unknown physical parameters; positivity is enforced by lower bounds.
    m.asymptote = pyo.Var(initialize=theta["asymptote"],
                         bounds=(1e-6, None), units=pyo.units.mg/pyo.units.L)
    m.rate_constant = pyo.Var(initialize=theta["rate_constant"],
                             bounds=(1e-6, None), units=1/pyo.units.day)
    # The input is fixed when fitting data and freed when designing a run.
    m.time_days = pyo.Var(initialize=data["time_days"],
                          bounds=design_bounds, units=pyo.units.day)
    # This variable is the prediction; the observed value is a suffix below.
    m.y = pyo.Var(initialize=data["y"], units=pyo.units.mg/pyo.units.L)
    # Estimation/design tools unfix the appropriate labeled variables.
    for v in (m.asymptote, m.rate_constant, m.time_days):
        v.fix()
    m.response = pyo.Constraint(expr=m.y == m.asymptote *
                                (1 - pyo.exp(-m.rate_constant * m.time_days)))
    return m

def label_bod_model(m, data, sigma):
    """Attach observations, unknowns, noise scale, and design inputs."""
    # Match each predicted output to its measured value.
    m.experiment_outputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
    m.experiment_outputs[m.y] = data["y"]
    # Tell ParmEst which parameters are shared across experiments.
    m.unknown_parameters = pyo.Suffix(direction=pyo.Suffix.LOCAL)
    for v in (m.asymptote, m.rate_constant):
        m.unknown_parameters[v] = pyo.value(v)
    # Store standard deviations, not variances, in output units (mg/L).
    m.measurement_error = pyo.Suffix(direction=pyo.Suffix.LOCAL)
    m.measurement_error[m.y] = sigma
    # Tell Pyomo.DoE which experimental input it may choose.
    m.experiment_inputs = pyo.Suffix(direction=pyo.Suffix.LOCAL)
    m.experiment_inputs[m.time_days] = data["time_days"]
    return m

Estimate the parameters¶

Create one experiment per row. The weighted sum of squared errors divides each residual by its measurement standard deviation before squaring.

ParmEst’s covariance interface requires more experiments than unknown parameters. With known noise variance, a full-rank information matrix can mathematically exist even when the two counts are equal; estimating noise from residuals additionally needs residual degrees of freedom.

Keep the objective normalization visible¶

SSE_weighted uses a per-experiment contribution $\tfrac12\sum_j(r_{ij}/\sigma_{ij})^2$ . ParmEst returns its average over experiments:

\bar\ell=\ell/N,\qquad \ell=N\bar\ell.

(9)

Multiplication by the positive constant $N$ leaves the fitted parameters unchanged. It does change curvature and likelihood-ratio statistics. This notebook forms nll_hat = len(data) * mean_nll before comparing a likelihood increase with a chi-squared cutoff. SSE instead uses unweighted squared residuals without the factor $1/2$ ; interpreting its value requires a separate noise-scale calculation.

def estimate_parameters(observations):
    """Fit the BOD model to a table of (time_days, y) rows.

    One `BODExperiment` is created per row. Returns the mean negative log
    likelihood, the fitted parameters, and the local covariance matrix -- everything
    needed both to report a fit and to seed the next experimental design.
    """
    experiments = [BODExperiment(row) for row in observations.to_dict("records")]
    estimator = Estimator(experiments, obj_function="SSE_weighted")

    weighted_sse, theta = estimator.theta_est()
    covariance = np.asarray(estimator.cov_est(method="finite_difference"))
    return weighted_sse, theta, covariance


weighted_sse, theta, covariance = estimate_parameters(data)

standard_deviation = np.sqrt(np.diag(covariance))
correlation = covariance / np.outer(standard_deviation, standard_deviation)

print(f"Mean half weighted SSE: {weighted_sse:.6f}")
print("Parameters:", theta)
print("Covariance:\n", covariance)
print("Correlation:\n", correlation)

assert np.isclose(theta["asymptote"], 19.1426, rtol=2e-4)
assert np.isclose(theta["rate_constant"], 0.53109, rtol=2e-4)

Mean half weighted SSE: 2.165856
Parameters: {'asymptote': 19.142575277697915, 'rate_constant': 0.5310913777423623}
Covariance:
 [[ 0.95876069 -0.06652725]
 [-0.06652725  0.00634735]]
Correlation:
 [[ 1.         -0.85280176]
 [-0.85280176  1.        ]]

# Archive both the measurements and the fitted curve for the handout.
fit_time = np.linspace(0, 8, 250)
results = {"time": data["time_days"].tolist(), "observed": data["y"].tolist(),
           "curve_time": fit_time.tolist(),
           "curve_y": (theta["asymptote"]*(1-np.exp(-theta["rate_constant"]*fit_time))).tolist(),
           "theta": dict(theta)}
helper.save_results("bod-fitted-response", results,
    notebook="notebooks/5-dev/BOD_Parmest_DoE.ipynb",
    source_tag="handout:l8-fit-curve", solver="Ipopt/MUMPS",
    description="Six BOD observations and the fitted nonlinear response",
        source_cells=[('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-observations'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-experiment-class'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-model'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-labels'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-estimate-parameters'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-fit'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-cov')])

[helper] wrote figures/results/bod-fitted-response.json

fig, ax = plt.subplots(figsize=(6.4, 3.5))
ax.plot(results["curve_time"], results["curve_y"], color="black", label="fitted response")
ax.scatter(results["time"], results["observed"], color="#D55E00", marker="o",
           label="observations", zorder=3)
ax.set_xlabel("Time (days)", fontsize=12)
ax.set_ylabel("Biochemical oxygen demand (mg/L)", fontsize=12)
ax.tick_params(labelsize=10)
ax.grid(alpha=0.25)
ax.legend(fontsize=10)
fig.tight_layout()
helper.save_figure(fig, "bod-fitted-response")
plt.show()

[helper] wrote media/figures/bod-fitted-response.png and .pdf

# MUMPS is the sparse linear solver bundled with the course Ipopt.
# This solver choice changes numerical implementation, not the likelihood.
experiments = [BODExperiment(row) for row in data.to_dict("records")]
pest = Estimator(experiments, obj_function="SSE_weighted",
                 solver_options={"linear_solver": "mumps"})
mean_nll, theta = pest.theta_est()
# ParmEst returns the average; the lecture uses the sum.
nll_hat = len(data) * mean_nll
print(f"theta = ({theta['asymptote']:.4f}, {theta['rate_constant']:.5f})")
print(f"Half weighted SSE (sum) = {nll_hat:.4f}")

theta = (19.1426, 0.53109)
Half weighted SSE (sum) = 12.9951

# Use expected information (Gauss-Newton) for the design demonstration.
cov_fd = pest.cov_est(method="finite_difference")
covariance = np.asarray(cov_fd)
# Other interfaces (see the validation note below before using them):
# cov_rh = pest.cov_est(method="reduced_hessian")
# cov_ad = pest.cov_est(method="automatic_differentiation_kaug")
print("Verified Gauss-Newton covariance (finite differences):")
print(cov_fd.round(5).to_string())

Verified Gauss-Newton covariance (finite differences):
               asymptote  rate_constant
asymptote        0.95876       -0.06653
rate_constant   -0.06653        0.00635

Sensitivities for explicit and implicit models¶

Stack the predictions in $\hat{\mathbf y}\in\mathbb R^{n_y}$ and use a row-oriented Jacobian:

\mathbf Q=\frac{d\hat{\mathbf y}}{d\boldsymbol\theta}\in\mathbb R^{n_y\times p}, \qquad Q_{aj}=\frac{d\hat y_a}{d\theta_j}.

(10)

For BOD, one observation contributes the sensitivity row

\mathbf q(t)=\begin{bmatrix}1-e^{-\theta_2t}&\theta_1t e^{-\theta_2t}\end{bmatrix}.

(11)

Both entries vanish at $t=0$ : such a measurement supplies no information about these two parameters under this response model. For large $t$ , sensitivity to the asymptote approaches one while sensitivity to the rate tends to zero.

For an implicit model, collect states and predictions in $\mathbf z$ and write $\mathbf F(\mathbf z(\boldsymbol\theta),\boldsymbol\varphi,\boldsymbol\theta)=\mathbf0$ . At a solution where $\mathbf F$ is continuously differentiable and the square matrix $\mathbf F_{\mathbf z}$ is nonsingular, the implicit-function theorem gives a locally unique differentiable solution branch. Differentiating at fixed design gives

\mathbf F_{\mathbf z}\mathbf Z+\mathbf F_{\boldsymbol\theta}=\mathbf0, \qquad \mathbf Z=\frac{d\mathbf z}{d\boldsymbol\theta}, \qquad \mathbf F_{\mathbf z}\mathbf Z=-\mathbf F_{\boldsymbol\theta}.

(12)

Solve this linear system for the sensitivity columns. If outputs are written separately as $\hat{\mathbf y}=\mathbf h(\mathbf x,\boldsymbol\varphi,\boldsymbol\theta)$ , include both paths of parameter dependence:

\mathbf Q=\mathbf h_{\boldsymbol\theta} +\mathbf h_{\mathbf x}\frac{d\mathbf x}{d\boldsymbol\theta}.

(13)

Automatic differentiation provides partial derivatives; differentiating through the model also requires the changing states. Forming an explicit matrix inverse is unnecessary. See Biegler (2010), §9.3.1, p. 261, Eqs. (9.21)--(9.22), and the erratum to the output-gradient formula on that page. A missing state/output dependence or a missing minus sign changes the sensitivity.

Total prediction sensitivity includes both direct and solved-state paths.

Differentiate the objective and obtain a Gauss--Newton step¶

Let $\mathbf W=\operatorname{blockdiag}(\boldsymbol\Sigma_1^{-1},\ldots,\boldsymbol\Sigma_N^{-1})$ be fixed. Because $\mathbf r=\mathbf y-\hat{\mathbf y}$ , its Jacobian is $-\mathbf Q$ . For $\ell=\tfrac12\mathbf r^T\mathbf W\mathbf r$ ,

\nabla\ell=-\mathbf Q^T\mathbf W\mathbf r, \qquad \nabla^2\ell=\mathbf Q^T\mathbf W\mathbf Q -\sum_{a=1}^{n_y}(\mathbf W\mathbf r)_a\nabla^2\hat y_a.

(14)

Gauss--Newton drops the residual-weighted output-curvature term. Equivalently, linearize the predictions at the current iterate and solve a weighted linear least-squares problem:

\hat{\mathbf y}(\boldsymbol\theta+\mathbf s) \approx\hat{\mathbf y}(\boldsymbol\theta)+\mathbf Q\mathbf s, \qquad \min_{\mathbf s}\frac12\|\mathbf W^{1/2}(\mathbf r-\mathbf Q\mathbf s)\|_2^2.

(15)

The normal equations and damped update are

(\mathbf Q^T\mathbf W\mathbf Q)\mathbf s=\mathbf Q^T\mathbf W\mathbf r, \qquad\boldsymbol\theta^+=\boldsymbol\theta+\alpha\mathbf s, \qquad0<\alpha\le1.

(16)

Full column rank makes the normal matrix positive definite. QR or SVD of the whitened Jacobian avoids explicitly forming an inverse and is preferable when conditioning is poor. The approximation is exact for an affine response. At a nonlinear noisy-data fit, stationarity means $\mathbf Q^T\mathbf W\mathbf r=\mathbf0$ , not that every residual is zero.

For the undamped iteration near a regular stationary solution, smoothness, full rank, and a sufficiently small omitted curvature term give local convergence. With zero residual and a Lipschitz Jacobian, the local rate is quadratic. A nonzero residual generally leaves a linear error term; sufficiently large residual-weighted curvature can prevent even local convergence. These statements require a close starting point and do not establish global convergence. Line searches and trust regions address failures of the full step. See Dennis & Schnabel, §10.2, Theorem 10.2.1 and Corollary 10.2.2, pp. 222--225; Seber & Wild, §2.1.3, pp. 25--27.

The independent BOD calculation below includes the omitted term in the observed Hessian and checks the Gauss--Newton covariance separately.

Gauss–Newton projects the observed response onto the local affine prediction model; the sketch uses an illustrative one-parameter response curve, not BOD data.

Check the covariance implementation¶

For this nonlinear model the inverse observed Hessian need not equal inverse expected information. The analytic calculation below validates the selected finite-difference method.

In the installed course build, cov_est(method="reduced_hessian") returned a matrix differing from the independently differentiated observed Hessian. Its call is shown above for API reference, but it is not used in this design demonstration. The third interface also requires a working k_aug installation. These software options are not interchangeable mathematical definitions.

# Independent checks of the mathematics and objective normalization.
t = data["time_days"].to_numpy(float)
a, k = theta["asymptote"], theta["rate_constant"]
e = np.exp(-k*t)
Q = np.column_stack([1-e, a*t*e])
r = data["y"].to_numpy(float) - a*(1-e)
M = Q.T @ Q  # sigma = 1
H = M.copy()
H[0, 1] -= np.sum(r*t*e)
H[1, 0] = H[0, 1]
H[1, 1] += np.sum(r*a*t*t*e)
assert np.allclose(cov_fd, np.linalg.inv(M), rtol=2e-5, atol=1e-7)
assert np.isclose(nll_hat, 0.5*(r@r), rtol=1e-10)
assert np.linalg.norm(Q.T @ r) < 1e-5
print("Independent inverse observed Hessian:")
print(np.linalg.inv(H))
from pyomo.util.check_units import assert_units_consistent
assert_units_consistent(experiments[0].get_labeled_model())

Independent inverse observed Hessian:
[[ 0.97048544 -0.06764595]
 [-0.06764595  0.00645409]]

import contextlib
import io
import logging

# Try all interfaces. Keep diagnostics available without flooding the handout.
covariance_methods = {"finite_difference": np.asarray(cov_fd)}
covariance_diagnostics = {}
for method in ("reduced_hessian", "automatic_differentiation_kaug"):
    diagnostic = io.StringIO()
    logging_level = logging.root.manager.disable
    try:
        logging.disable(logging.CRITICAL)
        with contextlib.redirect_stdout(diagnostic), contextlib.redirect_stderr(diagnostic):
            value = np.asarray(pest.cov_est(method=method), dtype=float)
        covariance_methods[method] = value
        print(method + ":")
        print(np.round(value, 6))
        # Compare each method with the curvature it is intended to compute.
        target = H if method == "reduced_hessian" else M
        target_name = "observed Hessian" if method == "reduced_hessian" else "information"
        print(f"Matches analytic inverse {target_name}:",
              np.allclose(value, np.linalg.inv(target), rtol=2e-4, atol=1e-6))
    except Exception as error:
        print(method + ": unavailable in this run (" + type(error).__name__ + ")")
        print(str(error).split("\n")[0][:180])
    finally:
        logging.disable(logging_level)
        covariance_diagnostics[method] = diagnostic.getvalue()
print("Analytic inverse observed Hessian (different target):")
print(np.round(np.linalg.inv(H), 6))

reduced_hessian:
[[ 0.947415 -0.065445]
 [-0.065445  0.006244]]
Matches analytic inverse observed Hessian: False

automatic_differentiation_kaug: unavailable in this run (ApplicationError)
Solver (asl) did not exit normally
Analytic inverse observed Hessian (different target):
[[ 0.970485 -0.067646]
 [-0.067646  0.006454]]

Expected information and estimator covariance¶

Before collecting data, $\mathbf Y$ is random. At a fixed parameter and design, the Gaussian measurement model gives $\mathbb E_{\boldsymbol\theta}[\mathbf r(\boldsymbol\theta;\mathbf Y)]=\mathbf0$ . The Jacobian and output Hessians are deterministic at that parameter, so

\begin{aligned} \mathbf M(\boldsymbol\theta) &=\mathbb E_{\boldsymbol\theta}[\nabla^2\ell(\boldsymbol\theta;\mathbf Y)]\\ &=\mathbf Q^T\mathbf W\mathbf Q -\sum_a\mathbb E_{\boldsymbol\theta}[(\mathbf W\mathbf r)_a]\nabla^2\hat y_a\\ &=\mathbf Q^T\mathbf W\mathbf Q =\sum_i\mathbf Q_i^T\boldsymbol\Sigma_i^{-1}\mathbf Q_i. \end{aligned}

(17)

This equality is exact for known, parameter-independent Gaussian covariance. It identifies expected Fisher information with the Gauss--Newton matrix, not with every realized observed Hessian. Under regularity conditions, information is also the covariance of the likelihood score. The score has the opposite sign to $\nabla\ell$ because $\ell$ is a negative log likelihood. See Held & Sabanés Bové (2014), §4.1.1, p. 81 and §5.4.

An estimator is unbiased if $\mathbb E_{\boldsymbol\theta}[\hat{\boldsymbol\theta}]=\boldsymbol\theta$ for every parameter value. It is consistent if $\hat{\boldsymbol\theta}_N\xrightarrow{P}\boldsymbol\theta$ as the sample grows. These are different properties; neither alone guarantees the other.

For an unbiased vector estimator, the Cramér--Rao inequality is

\operatorname{Cov}_{\boldsymbol\theta}(\hat{\boldsymbol\theta}) \succeq\mathbf M(\boldsymbol\theta)^{-1}.

(18)

This requires a regular statistical model, differentiation/expectation interchange, finite covariance, and nonsingular information at an interior parameter. The matrix inequality compares variance in every fixed direction. It does not assert that a nonlinear MLE is unbiased at finite sample size. Under appropriate identifiability and large-sample regularity conditions, MLEs are asymptotically normal and efficient. We estimate their covariance by

\widehat{\mathbf V}=\mathbf M(\hat{\boldsymbol\theta})^{-1}.

(19)

This is a local approximation for six BOD observations. Full rank at a fitted point alone does not establish the large-sample assumptions or correct model specification. See Held & Sabanés Bové, §4.2.1, pp. 95--98 and §5.4.2, pp. 144--145; Seber & Wild, Eq. (5.2), p. 191.

What a confidence region covers¶

A confidence procedure produces a random set $\mathcal C_{1-\alpha}(\mathbf Y)$ . Its coverage at a fixed true parameter is

\Pr_{\boldsymbol\theta}\{\boldsymbol\theta\in\mathcal C_{1-\alpha}(\mathbf Y)\}.

(20)

The target is $1-\alpha$ over repeated data sets. The true parameter stays fixed while the estimates and sets change. An observed frequentist confidence region does not assign a posterior probability to the parameter.

The random residual-curvature term averages to zero at the true parameter.

The local fit maps measurement errors into parameter errors; its propagated covariance is inverse information for a full-rank unconstrained affine model.

From standardized errors to a confidence ellipsoid¶

If $Z_1,\ldots,Z_k$ are independent standard normal variables, then $S=\sum_jZ_j^2$ has distribution $\chi_k^2$ , mean $k$ , and variance $2k$ . Its support is nonnegative. The 95th percentile leaves a 5% upper tail; it is not a normal quantile. See Held & Sabanés Bové, Appendix A.5, p. 337.

For a zero-mean multivariate normal error $\mathbf e=\hat{\boldsymbol\theta}-\boldsymbol\theta$ with covariance $\mathbf V\succ0$ , write $\mathbf V=\mathbf U\operatorname{diag}(v_j)\mathbf U^T$ , where $\mathbf U$ is orthogonal and $v_j>0$ . Rotation and scaling give

\mathbf u=\operatorname{diag}(v_j^{-1/2})\mathbf U^T\mathbf e, \qquad\operatorname{Cov}(\mathbf u)=\mathbf I, \qquad \|\mathbf u\|_2^2=\mathbf e^T\mathbf V^{-1}\mathbf e\sim\chi_p^2.

(21)

For a jointly normal vector, uncorrelated components are independent. Uncorrelated nonnormal components need not be independent. Applying the normal approximation to the fitted parameters and replacing $\mathbf V$ by its estimate makes this chi-squared calibration approximate.

from scipy.stats import chi2

# Deterministic probability densities: no random simulation or fitting is needed.
chi_x = np.linspace(0.02, 12, 800)
results = {"x": chi_x.tolist(), "df": [1, 2, 5],
           "density": {str(k): chi2.pdf(chi_x, k).tolist() for k in (1, 2, 5)},
           "cutoffs": {str(k): float(chi2.ppf(0.95, k)) for k in (1, 2, 5)}}
assert np.isclose(results["cutoffs"]["2"], 5.991464547107979)
helper.save_results("bod-chi-squared", results,
    notebook="notebooks/5-dev/BOD_Parmest_DoE.ipynb",
    source_tag="handout:l8-chi-squared", solver="scipy.stats.chi2 (no solver)",
    description="Chi-squared densities and 95% quantiles")

[helper] wrote figures/results/bod-chi-squared.json

fig, axes = plt.subplots(1, 2, figsize=(7.2, 3.2))
x = np.asarray(results["x"])
for k, color, style in zip(results["df"], ["#0072B2", "#D55E00", "#009E73"], ["-", "--", "-."]):
    axes[0].plot(x, results["density"][str(k)], color=color, linestyle=style,
                 label=f"{k} degree of freedom" if k == 1 else f"{k} degrees of freedom")
axes[0].set(xlim=(0, 12), ylim=(0, 0.65), xlabel="Squared standardized distance", ylabel="Probability density")
axes[0].legend(fontsize=8)
axes[0].set_title("Chi-squared densities")
y = np.asarray(results["density"]["2"])
c = results["cutoffs"]["2"]
axes[1].plot(x, y, color="black")
axes[1].fill_between(x, y, where=x<=c, color="#0072B2", alpha=0.2, label="95% area")
axes[1].fill_between(x, y, where=x>=c, color="#D55E00", alpha=0.5, hatch="///", label="5% upper tail")
axes[1].axvline(c, color="black", linestyle="--")
axes[1].text(c+0.15, 0.24, f"cutoff = {c:.3f}", fontsize=9, rotation=90)
axes[1].set(xlim=(0, 12), ylim=(0, 0.55), xlabel="Squared standardized distance", title="95% cutoff: two degrees of freedom")
axes[1].legend(fontsize=8, loc="upper right")
for ax in axes:
    ax.set_xlabel("Squared standardized distance", fontsize=10)
    ax.set_ylabel(ax.get_ylabel(), fontsize=10)
    ax.set_title(ax.get_title(), fontsize=11)
    ax.tick_params(labelsize=9)
    ax.grid(alpha=0.2)
fig.tight_layout()
helper.save_figure(fig, "bod-chi-squared")
plt.show()

[helper] wrote media/figures/bod-chi-squared.png and .pdf

Ellipsoid axes and the role of noise scale¶

Invert the standardized-distance test to obtain the local region

\mathcal C_{1-\alpha}^{\mathrm W} =\{\boldsymbol\theta:(\boldsymbol\theta-\hat{\boldsymbol\theta})^T \widehat{\mathbf V}^{-1}(\boldsymbol\theta-\hat{\boldsymbol\theta}) \le c\},\qquad c=\chi^2_{p,1-\alpha}.

(22)

In covariance eigenvector coordinates, the boundary is $\sum_j w_j^2/v_j=c$ . Its semi-axis along eigenvector $j$ has length $\sqrt{cv_j}$ . If $\widehat{\mathbf V}=\mathbf M^{-1}$ , the same axis has length $\sqrt{c/\lambda_j}$ for the information eigenvalue paired with that eigenvector. Inversion reverses ascending eigenvalue order. A long axis identifies a weakly informed combination in the chosen parameter units.

The $p$ -dimensional volume is $b_p c^{p/2}\sqrt{\det\widehat{\mathbf V}}$ , where $b_p$ is the unit-ball volume. For $p=2$ , the area is $\pi c\sqrt{\det\widehat{\mathbf V}}$ . The axis ratio is affected by parameter units; rescaling axes in a plotting program does not change the statistical construction.

Known versus estimated variance. With a known common $\sigma^2$ , $\widehat{\mathbf V}=\sigma^2(\mathbf Q^T\mathbf Q)^{-1}$ and the cutoff is $\chi^2_{p,1-\alpha}$ . If variance is estimated as $s^2=\mathrm{SSE}_{\min}/(n-p)$ , the usual linearization region instead has

(\boldsymbol\theta-\hat{\boldsymbol\theta})^T\mathbf Q^T\mathbf Q (\boldsymbol\theta-\hat{\boldsymbol\theta}) \le p s^2F_{p,n-p;1-\alpha}.

(23)

This F calibration is exact for a full-rank Gaussian linear model and approximate for nonlinear regression. Bates & Watts, Eq. (2.13), p. 53 and Seber & Wild, Eq. (5.9), p. 194 use this estimated-variance version. Held & Sabanés Bové, §5.4.2, p. 145 gives the multivariate normal/Wald construction. The known-noise BOD ellipse below uses $p=2$ and $c=5.9914645$ ; its nominal coverage is not a small-sample guarantee.

Rotate the correlated parameter errors and divide by directional standard deviations before using a central chi-squared cutoff.

# The local ellipse closes Part I before the nonlinear comparisons in Part II.
local_ellipse_results = {
    "theta": dict(theta), "covariance": covariance.tolist(),
    "cutoff": float(chi2.ppf(0.95, 2)), "sigma": MEASUREMENT_STD,
}
helper.save_results("bod-local-ellipse", local_ellipse_results,
    notebook="notebooks/5-dev/BOD_Parmest_DoE.ipynb",
    source_tag="handout:l8-local-ellipse-data",
    source_cells=[("notebooks/5-dev/BOD_Parmest_DoE.ipynb", tag) for tag in
                  ("evidence:l8-observations", "handout:l8-bod-model",
                   "handout:l8-bod-labels", "evidence:l8-estimate-parameters",
                   "handout:l8-cov")],
    solver="Ipopt/MUMPS; finite-difference information",
    description="Known-noise local 95% BOD confidence ellipse for Part I")
results = local_ellipse_results

[helper] wrote figures/results/bod-local-ellipse.json

def plot_local_ellipse(results):
    """Draw the actual BOD covariance ellipse, with its two semi-axes."""
    center = np.array([results["theta"]["asymptote"], results["theta"]["rate_constant"]])
    values, vectors = np.linalg.eigh(np.asarray(results["covariance"]))
    lengths = np.sqrt(results["cutoff"] * values)
    angle = np.linspace(0, 2*np.pi, 361)
    ellipse = center[:, None] + vectors @ (lengths[:, None] * np.vstack([np.cos(angle), np.sin(angle)]))
    fig, ax = plt.subplots(figsize=(6.4, 4.2))
    ax.plot(*ellipse, color="black", label="Local 95% boundary")
    for j, (color, style) in enumerate([("#0072B2", "--"), ("#D55E00", ":")]):
        ends = center[:, None] + np.outer(vectors[:, j] * lengths[j], [-1, 1])
        ax.plot(*ends, color=color, linestyle=style, linewidth=2,
                label=f"Axis {j+1}")
    ax.plot(*center, "kx", markersize=9, markeredgewidth=2, label="Fitted parameters")
    ax.set(xlabel=r"$\theta_1$ (mg/L)", ylabel=r"$\theta_2$ (day$^{-1}$)")
    ax.legend(fontsize=9)
    fig.tight_layout()
    return fig

fig = plot_local_ellipse(results)
helper.save_figure(fig, "bod-local-ellipse")
plt.show()

[helper] wrote media/figures/bod-local-ellipse.png and .pdf

Part II: nonlinear uncertainty and experimental design¶

Likelihood-ratio regions beyond the ellipse¶

For testing all $p$ parameters jointly, compare a candidate with the unrestricted MLE:

\Lambda(\boldsymbol\theta_0)=\frac{L(\boldsymbol\theta_0)}{L(\hat{\boldsymbol\theta})}, \qquad D(\boldsymbol\theta_0)=-2\log\Lambda =2[\ell(\boldsymbol\theta_0)-\ell(\hat{\boldsymbol\theta})].

(24)

The likelihood-ratio region retains candidates with $D\le\chi^2_{p,1-\alpha}$ . It evaluates the nonlinear likelihood, rather than imposing an elliptical shape. The cutoff still relies on regularity and large-sample calibration. With known common $\sigma$ , the BOD statistic is

D(\boldsymbol\theta) =\frac{\mathrm{SSE}(\boldsymbol\theta)-\mathrm{SSE}(\hat{\boldsymbol\theta})}{\sigma^2}.

(25)

The later comparison plots use the same original fit, the same six observations, and $\sigma=1$ mg/L for both the nonlinear region and the Gauss--Newton ellipse. A joint test fixes both BOD parameters and uses 2 degrees of freedom. A profile test fixes one and reoptimizes the other, so it uses 1. Parameter boundaries, unidentified parameters, or small samples can invalidate the usual chi-squared approximation. See Held & Sabanés Bové, §5.4.4, pp. 146--147, Eq. (5.15).

Why the textbook’s BOD contours are much larger¶

Bates & Watts, Example BOD 4, p. 51 estimates $s^2=6.498$ on 4 residual degrees of freedom. Its BOD contours in §6.1.1, p. 202 use estimated variance and an F-based construction. They extend to arbitrarily large rate constants. Those contours do not use this notebook’s assumed known $\sigma=1$ .

The difference can be checked directly. As $\theta_2\to\infty$ , all positive-time responses approach $\theta_1$ . The best horizontal line is $\theta_1=\bar y$ , with limiting SSE $\sum_i(y_i-\bar y)^2\approx107.2$ , compared with the nonlinear minimum $\approx25.99$ . The known-unit-variance statistic tends to $\approx81.21$ , far above the joint 95% cutoff 5.991. The same limiting fit can lie inside the much broader estimated-variance, F-calibrated region. Boundedness is conditional on the error model, noise scale, confidence level, and parameter domain.

Curvature distinguishes bending of the response manifold away from its tangent plane (intrinsic curvature) from distortion of coordinates along that manifold (parameter-effects curvature). Reparametrization can reduce the latter, while intrinsic bending remains. Curvature diagnostics and nonlinear regions help assess the local approximation; a small-looking ellipse by itself does not validate it. See Bates & Watts, Chapter 7, pp. 232, 245--246. No curvature measure is computed by the cells below.

Bootstrap: simulate or resample, then refit¶

A bootstrap approximates the sampling distribution by repeatedly refitting artificial data sets. For the fixed BOD times and assumed known Gaussian noise, use the parametric bootstrap below. A residual bootstrap instead resamples centered residuals at fixed times; a case bootstrap resamples whole (time, response) experiments. These answer different sampling questions.

The 200 replicates are a classroom demonstration, not a precision claim about interval endpoints. Bootstrap intervals are approximate and can be unreliable with six observations. The assumed standard deviation 1 is illustrative; the historical residuals are larger than this noise assumption predicts.

Regression resampling reference: Wehrens, Putter & Buydens (2000), §7.2, pp. 47–48, DOI.

Use fixed-design parametric resampling when the planned times and assumed error law represent the repeated experiment. Case resampling instead redraws complete (time, response) pairs; some times disappear and others repeat. It is appropriate for sampling cases from a population, and need not represent a fixed experimental design. See Wehrens, Putter, and Buydens (2000), Section 7.2, pp. 47–48. The bootstrap cloud estimates sampling variation; it is not itself a calibrated joint confidence region or a Bayesian posterior.

The artificial data and refits¶

First fit the original observations. For bootstrap replicate $b$ , keep the planned times and generate independent errors from the assumed measurement model:

y_i^{*(b)}=\hat y(t_i;\hat{\boldsymbol\theta})+\varepsilon_i^{*(b)}, \quad\varepsilon_i^{*(b)}\sim\mathcal N(0,\sigma^2), \qquad \hat{\boldsymbol\theta}^{*(b)}\in\arg\min_{\boldsymbol\theta\in\Theta} \frac12\sum_i\left[\frac{y_i^{*(b)}-\hat y(t_i;\boldsymbol\theta)}{\sigma}\right]^2.

(26)

Use the empirical distribution of these refits to estimate sampling variation. A percentile interval for component $j$ is $[q_{\alpha/2}(\hat\theta_j^*),q_{1-\alpha/2}(\hat\theta_j^*)]$ . Such componentwise intervals are not a simultaneous joint region. More replicates reduce Monte Carlo error, not model misspecification or limitations of six original observations.

For residual resampling with exchangeable homoscedastic errors, center the fitted residuals, sample them with replacement, and add them to the fitted response at each fixed time. Heteroscedastic or dependent errors require an appropriate alternative; blindly resampling a common residual pool discards that structure. Case resampling instead redraws whole experiments, changing which times occur. See Wehrens, Putter & Buydens (2000), §7.2, pp. 47--48. Inspect numerical failures and extreme fits before reporting intervals; do not silently drop them.

Fix the original fit during simulation, then refit every artificial data set.

# Parametric bootstrap: keep times fixed; simulate noise and refit.
rng = np.random.default_rng(804)
bootstrap_theta = []
for b in range(200):
    sample = data.copy()
    sample["y"] = bod_response(data["time_days"], theta) + rng.normal(0, 1, len(data))
    exps = [BODExperiment(row) for row in sample.to_dict("records")]
    fit = Estimator(exps, obj_function="SSE_weighted",
                    solver_options={"linear_solver": "mumps"})
    _, theta_b = fit.theta_est()
    bootstrap_theta.append(dict(theta_b))
bootstrap_theta = pd.DataFrame(bootstrap_theta)
print("Illustrative 95% percentile intervals (200 replicates):")
print(bootstrap_theta.quantile([0.025, 0.975]).round(4).to_string())

Illustrative 95% percentile intervals (200 replicates):
       asymptote  rate_constant
0.025     17.663         0.3958
0.975     21.479         0.6916

Built-in case bootstrap¶

theta_est_bootstrap resamples Experiment objects. It does not implement the fixed-time parametric algorithm above. ParmEst also filters insufficiently diverse resamples; with six observations this distinction matters.

# Resample whole (time, response) experiments; retain sample indices.
case_theta = pest.theta_est_bootstrap(100, seed=804, return_samples=True)
# Extreme values flag fits to inspect, not intervals to trust.
print("Case fits with rate > 10:", int((case_theta["rate_constant"] > 10).sum()))

Case fits with rate > 10: 8

Some case resamples produce extremely large rate estimates. The response becomes nearly saturated at all selected times. This can reflect weak information or a solver stopping at a poor stationary point. Independent checks found lower objective values for several returned fits. A solver reporting convergence does not make a fit reliable. Inspect the fit and try other initial values before constructing intervals; do not silently delete these samples. The fixed-design parametric bootstrap above is a different sampling experiment.

# Diagnose case fits without deleting or replacing samples.
# At fixed rate the least-squares asymptote has a closed-form solution.
poor_fits = 0
for _, fit in case_theta.iterrows():
    sample = data.iloc[fit["samples"]]
    ts, ys = sample["time_days"].to_numpy(), sample["y"].to_numpy()
    returned_sse = np.sum((ys - fit["asymptote"] *
                          (1-np.exp(-fit["rate_constant"]*ts)))**2)
    best_trial = returned_sse
    for trial_rate in (0.1, 0.2, 0.3, 0.5, 0.75, 1.0, 2.0):
        basis = 1-np.exp(-trial_rate*ts)
        trial_asymptote = (basis @ ys)/(basis @ basis)
        best_trial = min(best_trial, np.sum((ys-trial_asymptote*basis)**2))
    poor_fits += best_trial < returned_sse - 1e-3
print("Case fits with a demonstrably better feasible point:", poor_fits)

Case fits with a demonstrably better feasible point: 2

Profile likelihood: fix one parameter, refit the others¶

For each trial value $c$ of parameter $j$ , solve

\ell_{p,j}(c)=\min\{\ell(\boldsymbol\theta):\boldsymbol\theta\in\Theta,\ \theta_j=c\}, \qquad D_j(c)=2[\ell_{p,j}(c)-\ell(\hat{\boldsymbol\theta})].

(27)

A conditional slice holds the other parameters fixed. A profile reoptimizes them, so the profile objective cannot exceed the corresponding feasible conditional slice. This distinction captures parameter tradeoffs. A pointwise approximate 95% interval keeps $D_j(c)\le\chi^2_{1,0.95}=3.841$ ; it tests one restriction, not two.

Bates & Watts’ profile sum of squares $\widetilde S_j(c)$ is not our negative log likelihood. With known common $\sigma$ ,

\ell_{p,j}(c)=\frac{\widetilde S_j(c)}{2\sigma^2}, \qquad D_j(c)=\frac{\widetilde S_j(c)-S(\hat{\boldsymbol\theta})}{\sigma^2}.

(28)

Their signed profile $t$ function uses estimated residual scale $s$ :

\tau_j(c)=\operatorname{sign}(c-\hat\theta_j) \frac{\sqrt{\widetilde S_j(c)-S(\hat{\boldsymbol\theta})}}{s}.

(29)

They compare $|\tau_j|$ with a Student $t$ cutoff on $n-p$ degrees of freedom. For nonlinear regression that calibration is nominal, not an exact-coverage theorem. A straight profile $t$ plot indicates linear behavior in those coordinates; a profile trace plots the reoptimized nuisance parameters against $c$ . See Bates & Watts, §6.1.2, Eqs. (6.3)--(6.4), pp. 205--207.

The stable Pyomo model below uses summed half SSE, the same BOD response, known unit noise, and the positive parameter domain. It supplies all plotted profile results. Extend a grid that lacks a threshold crossing, inspect parameter bounds and solve status, and check alternative starts if a constrained fit finds an objective below the supposed unrestricted minimum. A truncated grid is not evidence of an unbounded interval.

Reference: Raue et al. (2009), pp. 1924--1926.

If the common noise variance is unknown¶

For $n$ independent scalar Gaussian errors, keep the variance-dependent term:

\ell(\boldsymbol\theta,\sigma^2)=\frac n2\log(\sigma^2)+\frac{S(\boldsymbol\theta)}{2\sigma^2}+C.

(30)

At fixed $\boldsymbol\theta$ with $S(\boldsymbol\theta)>0$ , minimization over variance gives $\hat\sigma^2(\boldsymbol\theta)=S(\boldsymbol\theta)/n$ , the ML variance estimate. Substitution gives $\ell_p(\boldsymbol\theta)=\tfrac n2\log[S(\boldsymbol\theta)/n]+n/2+C$ . Profiling the other response parameters as well yields

D_j(c)=n\log\!\left[\frac{\widetilde S_j(c)}{S(\hat{\boldsymbol\theta})}\right].

(31)

This is not $[\widetilde S_j(c)-S(\hat{\boldsymbol\theta})]/s^2$ . The ML variance estimate uses divisor $n$ , while the conventional residual estimate uses $n-p$ . A perfect fit $S=0$ makes the likelihood unbounded as variance tends to zero. The stable loop below assumes known unit variance; changing that assumption requires changing the statistic and its interpretation. See Held & Sabanés Bové, §5.3, pp. 130--131 for profiling a Gaussian variance nuisance parameter.

At fixed parameter 1, a profile refits parameter 2 along the vertical line. A slice keeps parameter 2 fixed. The homothetic ellipses are illustrative quadratic level sets, not the actual BOD likelihood.

def profile_model(observations):
    """BOD negative log likelihood, with shared parameters."""
    m = pyo.ConcreteModel()
    m.asymptote = pyo.Var(initialize=theta["asymptote"], bounds=(1e-6, None))
    m.rate_constant = pyo.Var(initialize=theta["rate_constant"], bounds=(1e-6, None))
    m.nll = pyo.Objective(expr=0.5 * sum(
        (row.y - m.asymptote * (1 - pyo.exp(-m.rate_constant * row.time_days)))**2
        for row in observations.itertuples()))
    return m


def profile_parameter(name, grid):
    """Fix one parameter and re-optimize the other at each grid point."""
    m = profile_model(data)
    solver = pyo.SolverFactory("ipopt")
    solver.options["linear_solver"] = "mumps"
    rows = []
    for value in grid:
        getattr(m, name).fix(float(value))
        pyo.assert_optimal_termination(solver.solve(m))
        rows.append((float(value), 2 * (pyo.value(m.nll) - nll_hat)))
    return np.asarray(rows)


profile_grids = {"asymptote": np.linspace(17, 22, 101),
                 "rate_constant": np.linspace(0.35, 0.8, 101)}
profiles = {name: profile_parameter(name, grid)
            for name, grid in profile_grids.items()}

Proposed ParmEst API (not executed)¶

Pyomo PR #3875 is open and draft as checked October 17, 2026, at head d86863725c676f71a4105c1f72060baa150b4e65. The installed course Pyomo 6.10.1 has no Estimator.profile_likelihood method. This preview requires that development branch and its dependencies and stays disabled here.

At this exact PR head, the method returns (profiles, intervals, theta_hat, obj_hat). For unregularized SSE_weighted, its lr_stat already equals 2 * len(exp_list) * (obj - obj_hat), converting the average objective to the summed likelihood scale. Do not multiply the returned statistic by the experiment count again. The SSE option instead profiles an unknown common variance and uses n_residuals * log(obj / obj_hat). These are different error-scale assumptions. See the pinned implementation.

The figure below is from the explicit stable-Pyomo implementation. The development preview has not been executed.

# DRAFT API: Pyomo PR #3875, head d868637 (not in Pyomo 6.10.1).
# Syntax preview only; results below use the stable Pyomo loop above.
RUN_DRAFT_API = False
if RUN_DRAFT_API:
    draft_profiles, draft_intervals, draft_theta_hat, draft_obj_hat = pest.profile_likelihood(
        profiled_theta=["asymptote", "rate_constant"],
        grid=profile_grids, theta_hat=dict(theta),
        obj_hat=mean_nll, warmstart="neighbor")
    # At this pinned head, lr_stat already includes the experiment count.
    # No additional normalization is applied.

from scipy.stats import chi2
cutoff = chi2.ppf(0.95, df=1)
for name, curve in profiles.items():
    # Interpolate the two threshold crossings; fail if the grid misses either.
    k = np.argmin(curve[:, 1])
    assert curve[0, 1] > cutoff and curve[-1, 1] > cutoff
    left = np.interp(cutoff, curve[:k+1, 1][::-1], curve[:k+1, 0][::-1])
    right = np.interp(cutoff, curve[k:, 1], curve[k:, 0])
    print(f"{name}: profile 95% interval [{left:.4f}, {right:.4f}]")
    assert curve[k, 1] >= -1e-5
print(f"Likelihood-ratio cutoff (1 parameter): {cutoff:.4f}")

asymptote: profile 95% interval [17.4284, 21.4132]
rate_constant: profile 95% interval [0.3915, 0.7133]
Likelihood-ratio cutoff (1 parameter): 3.8415

Identifiability depends on the model, designs, and data¶

Global structural identifiability over allowed designs means

[\mathbf f(\boldsymbol\varphi,\boldsymbol\theta) =\mathbf f(\boldsymbol\varphi,\widetilde{\boldsymbol\theta}) \text{ for every allowed }\boldsymbol\varphi] \ \Longrightarrow\ \boldsymbol\theta=\widetilde{\boldsymbol\theta}.

(32)

It concerns ideal predictions and the parameter domain. Practical identifiability concerns whether finite noisy observations constrain the parameters enough for the intended inference. A flat, one-sided, or unbounded profile can reveal weak practical information even when a local covariance can be computed. See Raue et al. (2009), §2.3, pp. 1924--1925.

At a fixed local fit, $\mathbf Q\mathbf v=\mathbf0$ for a nonzero $\mathbf v$ implies $\mathbf M\mathbf v=\mathbf0$ . This direction has no first-order effect on the available predictions. Full local rank does not rule out a different distant parameter vector giving identical predictions. Conversely, singular first-order information at an exceptional point need not prove global nonidentifiability; higher-order changes can matter.

The known- $\sigma=1$ profile intervals below have two threshold crossings. Bates & Watts’ BOD plots, pp. 206--207, use estimated noise and Student $t$ calibration: their 95% intervals are much wider, and their 99% intervals have no finite right endpoint. These results concern different inferential assumptions and levels.

results = {"bootstrap": bootstrap_theta.to_dict("list"),
           "profiles": {k: v.tolist() for k, v in profiles.items()},
           "cutoff": float(cutoff), "theta": dict(theta),
           "covariance": np.asarray(cov_fd).tolist()}
helper.save_results("bod-uncertainty", results,
    notebook="notebooks/5-dev/BOD_Parmest_DoE.ipynb",
    source_tag="handout:l8-profile-model", solver="Ipopt/MUMPS",
    description="BOD parametric bootstrap and explicit Pyomo likelihood profiles",
        source_cells=[('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-observations'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-experiment-class'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-model'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-labels'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-estimate-parameters'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-fit'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-cov'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bootstrap'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-profile-solver')])

[helper] wrote figures/results/bod-uncertainty.json

def plot_bod_uncertainty(results):
    """Compare a bootstrap cloud and likelihood profiles."""
    fig, axes = plt.subplots(1, 3, figsize=(11, 3.2))
    boot = results["bootstrap"]
    axes[0].scatter(boot["asymptote"], boot["rate_constant"],
                    s=10, alpha=0.5, color="#0072B2", marker="o")
    axes[0].plot(results["theta"]["asymptote"], results["theta"]["rate_constant"],
                 "kx", ms=8)
    axes[0].set(xlabel=r"$\theta_1$", ylabel=r"$\theta_2$", title="Parametric bootstrap")
    for ax, (name, values) in zip(axes[1:], results["profiles"].items()):
        curve = np.asarray(values)
        ax.plot(curve[:, 0], curve[:, 1], color="#0072B2", label="profile")
        ax.axhline(results["cutoff"], color="black", ls="--", label="95% cutoff")
        ax.set(xlabel=name, ylabel=r"$2(\ell_{p}-\ell_{\min})$", ylim=(0, 8))
        ax.legend(fontsize=8)
    fig.tight_layout()
    return fig


fig = plot_bod_uncertainty(results)
helper.save_figure(fig, "bod-uncertainty")

[helper] wrote media/figures/bod-uncertainty.png and .pdf

profile_results = results  # Preserve the one-parameter profile archive.
from scipy.stats import chi2
# Joint likelihood region: evaluate the nonlinear objective on a 2D grid.
asymptote_grid = np.linspace(15.5, 23.5, 241)
rate_grid = np.linspace(0.25, 0.95, 241)
AA, KK = np.meshgrid(asymptote_grid, rate_grid)
predicted = AA[..., None] * (1 - np.exp(-KK[..., None] * data["time_days"].to_numpy()))
joint_lr = np.sum((data["y"].to_numpy() - predicted)**2, axis=-1) - 2*nll_hat
joint_cutoff = float(chi2.ppf(0.95, 2))
# Boundary must be contained in the grid, rather than silently clipped.
assert min(joint_lr[0].min(), joint_lr[-1].min(),
           joint_lr[:, 0].min(), joint_lr[:, -1].min()) > joint_cutoff
region_results = {"asymptote_grid": asymptote_grid.tolist(),
                  "rate_grid": rate_grid.tolist(), "joint_lr": joint_lr.tolist(),
                  "cutoff": joint_cutoff, "theta": dict(theta),
                  "covariance": covariance.tolist(),
                  "bootstrap": bootstrap_theta.to_dict("list")}
for name in ("bod-confidence-regions", "bod-bootstrap-regions"):
    helper.save_results(name, region_results,
        notebook="notebooks/5-dev/BOD_Parmest_DoE.ipynb",
        source_tag="handout:l8-regions-data", solver="Ipopt/MUMPS + direct likelihood evaluation",
        description="Joint nonlinear likelihood region and GN ellipse for the same BOD fit",
        source_cells=[('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-observations'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-experiment-class'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-model'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-labels'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-estimate-parameters'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-fit'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-cov'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bootstrap'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-profile-solver')])

results = region_results

[helper] wrote figures/results/bod-confidence-regions.json

[helper] wrote figures/results/bod-bootstrap-regions.json

def plot_region_comparison(results):
    """Compare nonlinear likelihood and local Gaussian regions."""
    fig, ax = plt.subplots(figsize=(7.2, 4.8))
    center = np.array([results["theta"]["asymptote"], results["theta"]["rate_constant"]])
    values, vectors = np.linalg.eigh(np.asarray(results["covariance"]))
    angle = np.linspace(0, 2*np.pi, 361)
    ellipse = center[:, None] + vectors @ (np.sqrt(results["cutoff"]*values)[:, None]
               * np.vstack([np.cos(angle), np.sin(angle)]))
    ax.contour(results["asymptote_grid"], results["rate_grid"], results["joint_lr"],
               levels=[results["cutoff"]], colors=["#0072B2"], linewidths=2.5)
    ax.plot([], [], color="#0072B2", ls="-", label="Nonlinear joint LR boundary")
    ax.plot(*ellipse, color="#D55E00", ls="--", lw=2.5, label="Gauss-Newton ellipse")
    ax.plot(*center, "kx", ms=9, mew=2, label="Original estimate")
    ax.set(xlabel=r"$\theta_1$ (mg/L)", ylabel=r"$\theta_2$ (day$^{-1}$)")
    ax.legend(fontsize=10, loc="upper right")
    fig.tight_layout()
    return fig

fig = plot_region_comparison(results)
helper.save_figure(fig, "bod-confidence-regions")

[helper] wrote media/figures/bod-confidence-regions.png and .pdf

def plot_region_comparison(results):
    """Compare nonlinear likelihood and local Gaussian regions."""
    fig, ax = plt.subplots(figsize=(7.2, 4.8))
    center = np.array([results["theta"]["asymptote"], results["theta"]["rate_constant"]])
    values, vectors = np.linalg.eigh(np.asarray(results["covariance"]))
    angle = np.linspace(0, 2*np.pi, 361)
    ellipse = center[:, None] + vectors @ (np.sqrt(results["cutoff"]*values)[:, None]
               * np.vstack([np.cos(angle), np.sin(angle)]))
    ax.contour(results["asymptote_grid"], results["rate_grid"], results["joint_lr"],
               levels=[results["cutoff"]], colors=["#0072B2"], linewidths=2.5)
    ax.plot([], [], color="#0072B2", ls="-", label="Nonlinear joint LR boundary")
    ax.plot(*ellipse, color="#D55E00", ls="--", lw=2.5, label="Gauss-Newton ellipse")
    ax.scatter(results["bootstrap"]["asymptote"], results["bootstrap"]["rate_constant"],
               s=15, alpha=0.45, color="0.35", marker="o", label="200 parametric refits")
    ax.plot(*center, "kx", ms=9, mew=2, label="Original estimate")
    ax.set(xlabel=r"$\theta_1$ (mg/L)", ylabel=r"$\theta_2$ (day$^{-1}$)")
    ax.legend(fontsize=10, loc="upper right")
    fig.tight_layout()
    return fig

fig = plot_region_comparison(results)
helper.save_figure(fig, "bod-bootstrap-regions")

[helper] wrote media/figures/bod-bootstrap-regions.png and .pdf

results = profile_results
helper.save_results("bod-profile-results", results,
    notebook="notebooks/5-dev/BOD_Parmest_DoE.ipynb",
    source_tag="handout:l8-profile-model", solver="Ipopt/MUMPS",
    description="One-parameter profile likelihood results",
        source_cells=[('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-observations'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-experiment-class'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-model'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-labels'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-estimate-parameters'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-fit'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-cov'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bootstrap'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-profile-solver')])

[helper] wrote figures/results/bod-profile-results.json

def plot_profiles(results):
    """Display each one-parameter likelihood profile at its own cutoff."""
    fig, axes = plt.subplots(1, 2, figsize=(9, 3.5))
    labels = [r"$\theta_1$ (mg/L)", r"$\theta_2$ (day$^{-1}$)"]
    for ax, (name, values), label in zip(axes, results["profiles"].items(), labels):
        curve = np.asarray(values)
        ax.plot(curve[:, 0], curve[:, 1], color="#0072B2", ls="-", label="Profile LR")
        ax.axhline(results["cutoff"], color="black", ls="--", label="95% cutoff (1 df)")
        ax.set(xlabel=label, ylabel=r"$2(\ell_p-\ell_{min})$", ylim=(0, 8))
        ax.legend(fontsize=9)
    fig.tight_layout()
    return fig

fig = plot_profiles(results)
helper.save_figure(fig, "bod-profile-results")

[helper] wrote media/figures/bod-profile-results.png and .pdf

From covariance to experimental design¶

At the fitted parameter and the stated error scale, the analytically checked finite-difference covariance gives the existing-data information $\mathbf M_{\mathrm{old}}=\widehat{\mathbf V}^{-1}$ . The software name prior_FIM refers to information from previous experiments here; it is not a Bayesian prior distribution.

For one independent new BOD measurement, use the row sensitivity

\mathbf q(t)=\begin{bmatrix}1-e^{-\hat\theta_2t}&\hat\theta_1t e^{-\hat\theta_2t}\end{bmatrix}, \qquad \mathbf M_{\mathrm{total}}(t)=\mathbf M_{\mathrm{old}}+\frac{\mathbf q(t)^T\mathbf q(t)}{\sigma^2}.

(33)

The column vector in the updated_fim implementation below is the transpose of this row; np.outer produces the same $2\times2$ contribution. A scalar observation adds at most rank-one information. Its direction matters as well as the sizes of its sensitivity components.

A local design minimizes a chosen scalar criterion $J$ at a fixed nominal parameter:

\min_{\boldsymbol\varphi\in\Phi} J[\mathbf M_{\mathrm{old}}(\hat{\boldsymbol\theta})+ \mathbf M_{\mathrm{new}}(\boldsymbol\varphi,\hat{\boldsymbol\theta})].

(34)

For the D-optimal BOD design, maximize $\log\det\mathbf M_{\mathrm{total}}(t)$ on $0.5\le t\le5$ days. The future measured value is unnecessary for choosing this local design; it is needed for the subsequent refit. After measurement, update both the fitted parameters and information. See Bates & Watts, §3.14.4, p. 127; Seber & Wild, §5.13.1, pp. 250--251; Laky et al. (2026), §2.1--2.2.

Independent measurement information adds at the same nominal parameter, coordinate scale and noise model.

t_design = np.linspace(0.5, 5, 300)
q_asymptote = 1 - np.exp(-theta["rate_constant"] * t_design)
q_rate = theta["asymptote"] * t_design * np.exp(-theta["rate_constant"] * t_design)


sensitivity_results = {"time": t_design.tolist(),
    "asymptote": q_asymptote.tolist(), "rate_constant": q_rate.tolist()}
helper.save_results("bod-sensitivities", sensitivity_results,
    notebook="notebooks/5-dev/BOD_Parmest_DoE.ipynb",
    source_tag="handout:l8-sensitivities",
    source_cells=[("notebooks/5-dev/BOD_Parmest_DoE.ipynb", "evidence:l8-estimate-parameters"),
                  ("notebooks/5-dev/BOD_Parmest_DoE.ipynb", "evidence:l8-observations")],
    solver="analytic BOD derivatives at the fitted parameters",
    description="BOD prediction sensitivities over the allowed next-measurement window")
results = sensitivity_results

[helper] wrote figures/results/bod-sensitivities.json

t_design = np.asarray(results["time"])
q_asymptote = np.asarray(results["asymptote"])
q_rate = np.asarray(results["rate_constant"])

fig, axes = plt.subplots(1, 2, figsize=(10, 3.5))
axes[0].plot(t_design, q_asymptote, color="black")
axes[0].set_ylabel(r"$\partial y/\partial\theta_1$")
axes[1].plot(t_design, q_rate, color="#0072B2", linestyle="--")
axes[1].set_ylabel(r"$\partial y/\partial\theta_2$")
for ax in axes:
    ax.set_xlabel("Measurement time (days)")
    ax.grid(True)
plt.tight_layout()
helper.save_figure(fig, "bod-sensitivities")
plt.show()

[helper] wrote media/figures/bod-sensitivities.png and .pdf

Optimize the next measurement¶

Pyomo 6.10 configures the objective and solver behavior in the DesignOfExperiments constructor. run_doe() performs the optimization and stores the summary in doe.results.

The bounds in BODExperiment restrict the new measurement to 0.5--5 days.

prior_fim = np.linalg.inv(covariance)
design_experiment = BODExperiment(
    {"time_days": 2.0, "y": 10.0},
    measurement_standard_deviation=1.0,
    theta=dict(theta),
    design_bounds=(0.5, 5.0),
)

# Use the MUMPS linear solver available in the course environment.
ipopt = pyo.SolverFactory("ipopt")
ipopt.options["linear_solver"] = "mumps"

doe = DesignOfExperiments(
    experiment=design_experiment,
    fd_formula="central",
    step=1e-3,
    objective_option=ObjectiveLib.determinant,
    prior_FIM=prior_fim,
    improve_cholesky_roundoff_error=True,
    solver=ipopt,
    tee=False,
)
doe.run_doe(model=doe.create_doe_model())

optimal_time = doe.results["Experiment Design"][0]
print(f"D-optimal next measurement: {optimal_time:.6f} days")

D-optimal next measurement: 1.499740 days

Verify with an independent calculation¶

For this two-parameter model, a dense grid is cheap. It provides an independent check on the software formulation and local optimizer.

def updated_fim(t):
    sensitivity = np.array(
        [
            1 - np.exp(-theta["rate_constant"] * t),
            theta["asymptote"] * t * np.exp(-theta["rate_constant"] * t),
        ]
    )
    return prior_fim + np.outer(sensitivity, sensitivity)


grid = np.linspace(0.5, 5, 2001)
grid_logdet = np.array([np.linalg.slogdet(updated_fim(t))[1] for t in grid])
grid_time = grid[np.argmax(grid_logdet)]

print(f"Dense-grid optimum: {grid_time:.6f} days")
print(f"Prior log determinant:   {np.linalg.slogdet(prior_fim)[1]:.6f}")
print(f"Updated log determinant: {grid_logdet.max():.6f}")
assert abs(optimal_time - grid_time) < 0.02

plt.plot(grid, grid_logdet, color="black")
plt.axvline(optimal_time, color="#D55E00", linestyle="--", label="Pyomo.DoE optimum")
plt.xlabel("Measurement time (days)")
plt.ylabel("Updated log determinant")
plt.grid(True)
plt.legend()
plt.show()

Dense-grid optimum: 1.499000 days
Prior log determinant:   6.401107
Updated log determinant: 6.742537

Information eigenvalues and design criteria¶

For $\mathbf M=\mathbf U\operatorname{diag}(\lambda_j)\mathbf U^T\succ0$ , orthogonality gives

\operatorname{tr}\mathbf M=\sum_j\lambda_j,\quad \det\mathbf M=\prod_j\lambda_j,\quad \operatorname{tr}\mathbf M^{-1}=\sum_j\lambda_j^{-1},\quad \kappa_2(\mathbf M)=\lambda_{\max}/\lambda_{\min}.

(35)

Criterion	Objective	What it controls in the local approximation	Pyomo.DoE option
A	Minimize $\sum_j1/\lambda_j$	Sum of parameter variances in the chosen units	`trace`
D	Maximize $\prod_j\lambda_j$	Minimize ellipsoid volume at a common cutoff	`determinant`
E	Maximize $\lambda_{\min}$	Minimize the largest covariance eigenvalue and longest semi-axis	`minimum_eigenvalue`
ME	Minimize $\lambda_{\max}/\lambda_{\min}$	Balance axes; it does not control overall size	`condition_number`
pseudo-A	Maximize $\sum_j\lambda_j$	Increase total information; not A-optimality	`pseudo_trace`

The ellipse’s longest-to-shortest semi-axis ratio is $\sqrt{\kappa_2(\mathbf M)}$ , not $\kappa_2(\mathbf M)$ . A-optimality minimizes trace of covariance, not the volume of an enclosing coordinate box. Such a box depends on covariance diagonal entries and its volume involves their product. Pseudo-A is also called a trace or T criterion; a finite trace can coexist with a zero eigenvalue and an unidentifiable direction.

For information eigenvalues $(16,1)$ and $(4,3)$ , the A objectives are $17/16$ and $7/12$ ; D objectives are 16 and 12; E objectives are 1 and 3; ME objectives are 16 and $4/3$ ; pseudo-A objectives are 17 and 7. Thus A, E and ME prefer the second design, while D and pseudo-A prefer the first. There is no universal ordering of designs independent of the objective. See Laky et al. (2026), §2.2, Eqs. (12)--(18) and Table 1.

The five curves below use the same original observations, fitted parameters, units, and allowed time interval. Their grid optima need not coincide.

# Evaluate all five criteria on the same time grid and parameter scaling.
criteria_time = np.linspace(0.5, 5, 901)
criteria_information = np.asarray([updated_fim(t) for t in criteria_time])
eigenvalues = np.linalg.eigvalsh(criteria_information)
assert np.all(eigenvalues > 0)
criteria_values = {
    "A: minimize trace inverse": np.sum(1/eigenvalues, axis=1),
    "D: maximize determinant": np.prod(eigenvalues, axis=1),
    "E: maximize smallest eigenvalue": eigenvalues[:, 0],
    "ME: minimize condition number": eigenvalues[:, 1]/eigenvalues[:, 0],
    "Pseudo-A: maximize trace": np.sum(eigenvalues, axis=1),
}
criteria_results = {"time": criteria_time.tolist(),
    "criteria": {key: value.tolist() for key, value in criteria_values.items()}}
helper.save_results("bod-design-criteria", criteria_results,
    notebook="notebooks/5-dev/BOD_Parmest_DoE.ipynb",
    source_tag="handout:l8-criteria-data", solver="analytic sensitivity grid",
    description="All five BOD design criteria at common nominal parameters and units",
        source_cells=[('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-observations'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-experiment-class'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-model'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-labels'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-estimate-parameters'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-fit'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-cov'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-design'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-design-grid')])
for label, values in criteria_values.items():
    index = np.argmin(values) if "minimize" in label else np.argmax(values)
    print(f"{label}: grid optimum at {criteria_time[index]:.3f} days")

results = criteria_results

[helper] wrote figures/results/bod-design-criteria.json
A: minimize trace inverse: grid optimum at 5.000 days
D: maximize determinant: grid optimum at 1.500 days
E: maximize smallest eigenvalue: grid optimum at 5.000 days
ME: minimize condition number: grid optimum at 5.000 days
Pseudo-A: maximize trace: grid optimum at 1.885 days

def plot_all_criteria(results):
    """Use individual ordinate scales so all five objectives remain legible."""
    fig, axes = plt.subplots(3, 2, figsize=(10, 9))
    t = np.asarray(results["time"])
    for ax, (label, values) in zip(axes.flat, results["criteria"].items()):
        values = np.asarray(values)
        index = np.argmin(values) if "minimize" in label else np.argmax(values)
        ax.plot(t, values, color="#0072B2", ls="-")
        ax.plot(t[index], values[index], "ko", ms=6)
        ax.axvline(t[index], color="0.35", ls=":", lw=1.5)
        ax.set_title(label, fontsize=12)
        ax.set_xlabel("New measurement time (days)", fontsize=12)
        ax.tick_params(labelsize=11)
        ax.text(0.98, 0.93, f"grid optimum: {t[index]:.3f} d",
                transform=ax.transAxes, ha="right", va="top", fontsize=10,
                bbox={"facecolor": "white", "edgecolor": "none", "alpha": 0.85})
    axes.flat[-1].axis("off")
    axes.flat[-1].text(0.05, 0.65, "Same six existing observations\nSame nominal parameter estimate\nSame units and scaling\nOne proposed new observation",
                       transform=axes.flat[-1].transAxes, va="top", fontsize=12, linespacing=1.6)
    fig.tight_layout()
    return fig

fig = plot_all_criteria(results)
helper.save_figure(fig, "bod-design-criteria")

[helper] wrote media/figures/bod-design-criteria.png and .pdf

Parameter scaling and Cholesky implementation¶

Let $\boldsymbol{\theta}=\mathbf{D}\boldsymbol{\eta}$ , where $\mathbf{D}=\operatorname{diag}(s_1,\ldots,s_p)$ and $s_j>0$ . The chain rule gives $\mathbf{Q}_\eta=\mathbf{Q}_\theta \mathbf{D}$ and $\mathbf{M}_\eta=\mathbf{D} \mathbf{M}_\theta \mathbf{D}$ . Thus parameter scaling changes the units in which uncertainty is compared and can change A-, E-, ME-, and pseudo-A-optimal designs. D-optimal rankings are unchanged by a fixed nonsingular scaling because $\det(\mathbf{M}_\eta)=\det(\mathbf{D})^2\det(\mathbf{M}_\theta)$ .

Multiplying the entire information matrix by one positive constant does not change its condition number: $\kappa_2(c\mathbf{M})=c\lambda_{max}/(c\lambda_{min})=\kappa_2(\mathbf{M})$ . scale_constant_value controls a common scale; scale_nominal_param_value multiplies sensitivity columns by their nominal parameter values. Do not confuse these two scaling operations.

For positive-definite $\mathbf{M}$ , a Cholesky factorization $\mathbf{M}=\mathbf{LL}^T$ has $L_{jj}>0$ and

\det \mathbf{M}=(\det \mathbf{L})^2=\left(\prod_j L_{jj}\right)^2,\qquad \log\det \mathbf{M}=2\sum_j\log L_{jj}.

(36)

Pyomo.DoE introduces $\mathbf{L}$ as variables and $\mathbf{M}=\mathbf{LL}^T$ as constraints for its smooth log-determinant formulation. Eigenvalue criteria use external calculations and derivatives. See Laky et al. (2026), Section 2 and Table 1.

A nonuniform parameter scaling can improve or worsen the condition number; improvement is not guaranteed. Scale the existing-data information consistently with new-experiment information before adding the matrices.

Why adding information reduces predicted covariance¶

At a fixed nominal parameter, write the new information as $\mathbf A^T\mathbf A\succeq0$ and suppose $\mathbf M_0\succ0$ . Let $\mathbf B=\mathbf A\mathbf M_0^{-1/2}$ . Then

\begin{aligned} \mathbf M_1&=\mathbf M_0+\mathbf A^T\mathbf A =\mathbf M_0^{1/2}(\mathbf I+\mathbf B^T\mathbf B)\mathbf M_0^{1/2},\\ \mathbf M_1^{-1} &=\mathbf M_0^{-1/2}(\mathbf I+\mathbf B^T\mathbf B)^{-1}\mathbf M_0^{-1/2} \preceq\mathbf M_0^{-1}. \end{aligned}

(37)

The eigenvalues of $(\mathbf I+\mathbf B^T\mathbf B)^{-1}$ lie in $(0,1]$ , so predicted variance cannot increase in any fixed direction. Equality may remain in directions the new experiment does not inform. Recomputing sensitivities at a different refitted parameter changes the premise: an observed covariance need not match this fixed-parameter prediction.

For one scalar measurement, set $\mathbf u=\mathbf q^T/\sigma$ and use the rank-one formulas

\frac{\det(\mathbf M_0+\mathbf u\mathbf u^T)}{\det\mathbf M_0} =1+\mathbf u^T\mathbf M_0^{-1}\mathbf u, \qquad (\mathbf M_0+\mathbf u\mathbf u^T)^{-1} =\mathbf M_0^{-1} -\frac{\mathbf M_0^{-1}\mathbf u\mathbf u^T\mathbf M_0^{-1}} {1+\mathbf u^T\mathbf M_0^{-1}\mathbf u}.

(38)

The D-optimal next measurement therefore maximizes its sensitivity in the existing covariance metric, $\mathbf q\widehat{\mathbf V}\mathbf q^T/\sigma^2$ . Merely choosing the largest Euclidean sensitivity norm ignores information already collected. These identities apply to the local additive-information model used here; they are not guarantees about the realized refit.

Close the loop: take the designed measurement¶

The design step produced a time, not data. Closing the loop means going back to the experiment, measuring at $t^\star$ , appending the result, and re-estimating.

This is exactly where the digital experiment is needed: we evaluate the true response at $t^\star$ , add one noise draw, and hand the seven-point data set back to parmest.

The quantity to watch is $\det\mathbf{V}_\theta$ . Because the 95% confidence region is the ellipse

(\boldsymbol{\theta}-\hat{\boldsymbol{\theta}})^{\mathsf{T}}\mathbf{V}_\theta^{-1}(\boldsymbol{\theta}-\hat{\boldsymbol{\theta}})\le\chi^2_{2,0.95}=5.991,

(39)

its area is $\pi\,\chi^2_{2,0.95}\sqrt{\det\mathbf{V}_\theta}$ . A smaller determinant is a smaller ellipse, and D-optimal design maximizes $\det\mathbf{F}$ , which is the same objective seen from the other side.

# A separate random stream for the closed loop, so the new measurement does not
# depend on how many draws the calibration table above happened to consume.
loop_rng = np.random.default_rng(2026)

new_observation = run_digital_experiment(optimal_time, loop_rng)
print(f"Measured at t* = {optimal_time:.4f} days:  y = {new_observation:.4f}")
print(f"Noise-free truth at t*:        y = {bod_response(optimal_time, THETA_TRUE):.4f}")

augmented = pd.concat(
    [data, pd.DataFrame({"time_days": [optimal_time], "y": [new_observation]})],
    ignore_index=True,
)
augmented

Measured at t* = 1.4997 days:  y = 9.6257
Noise-free truth at t*:        y = 10.4188

CHI2_95_2DOF = 5.991  # 95% quantile of the chi-squared distribution, 2 degrees of freedom


def ellipse_area(cov):
    """Area of the 95% confidence ellipse implied by a 2x2 covariance matrix."""
    return np.pi * CHI2_95_2DOF * np.sqrt(np.linalg.det(cov))


def summarize(label, theta, cov):
    """One row of the before/after comparison table."""
    sd = np.sqrt(np.diag(cov))
    return {
        "design": label,
        "theta_1": theta["asymptote"],
        "theta_2": theta["rate_constant"],
        "sd(theta_1)": sd[0],
        "sd(theta_2)": sd[1],
        "det(V)": np.linalg.det(cov),
        "ellipse area": ellipse_area(cov),
    }


sse_after, theta_after, covariance_after = estimate_parameters(augmented)

print("Covariance before (6 observations):\n", covariance)
print("\nCovariance after  (7 observations, measured at t*):\n", covariance_after)

closed_loop = pd.DataFrame(
    [
        summarize("before: 6 observations", theta, covariance),
        summarize(f"after: + measurement at t* = {optimal_time:.3f}", theta_after, covariance_after),
    ]
)
closed_loop

Covariance before (6 observations):
 [[ 0.95876069 -0.06652725]
 [-0.06652725  0.00634735]]

Covariance after  (7 observations, measured at t*):
 [[ 0.98389822 -0.05729307]
 [-0.05729307  0.00447954]]

determinant_drop = 1 - np.linalg.det(covariance_after) / np.linalg.det(covariance)
area_drop = 1 - ellipse_area(covariance_after) / ellipse_area(covariance)

print(f"det(V) fell by {100 * determinant_drop:.1f}%")
print(f"95% ellipse area fell by {100 * area_drop:.1f}%")

# Sanity check: one well-placed measurement should measurably shrink the region.
assert np.linalg.det(covariance_after) < np.linalg.det(covariance)

det(V) fell by 32.2%
95% ellipse area fell by 17.7%

sequential_results = {"before_theta": dict(theta), "after_theta": dict(theta_after),
    "before_cov": covariance.tolist(), "after_cov": covariance_after.tolist(),
    "cutoff": float(chi2.ppf(0.95, 2)), "time": float(optimal_time)}
helper.save_results("bod-sequential-ellipses", sequential_results,
    notebook="notebooks/5-dev/BOD_Parmest_DoE.ipynb",
    source_tag="handout:l8-sequential-data", solver="Ipopt/MUMPS",
    description="Actual and centered uncertainty regions before and after seventh observation",
        source_cells=[('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-observations'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-experiment-class'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-model'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-bod-labels'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-estimate-parameters'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-fit'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-cov'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-design'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-digital-experiment'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'handout:l8-measure'), ('notebooks/5-dev/BOD_Parmest_DoE.ipynb', 'evidence:l8-refit')])

results = sequential_results

[helper] wrote figures/results/bod-sequential-ellipses.json

def plot_sequential_regions(results):
    """Separate the actual estimate shift from the change in ellipse shape."""
    fig, axes = plt.subplots(1, 2, figsize=(10, 4.2))
    angles = np.linspace(0, 2*np.pi, 361)
    circle = np.vstack([np.cos(angles), np.sin(angles)])
    before = np.array([results["before_theta"]["asymptote"], results["before_theta"]["rate_constant"]])
    after = np.array([results["after_theta"]["asymptote"], results["after_theta"]["rate_constant"]])
    for ax, centered in zip(axes, [False, True]):
        for center, cov, style, color, label in [
            (before, results["before_cov"], "--", "#D55E00", "6 observations"),
            (after, results["after_cov"], "-", "#0072B2", "7 observations")]:
            values, vectors = np.linalg.eigh(np.asarray(cov))
            shown_center = before if centered else center
            points = shown_center[:, None] + vectors @ (np.sqrt(results["cutoff"]*values)[:, None]*circle)
            ax.plot(*points, ls=style, color=color, label=label)
            ax.plot(*shown_center, marker="x" if label.startswith("6") else "+", color=color, ms=8)
        ax.set(xlabel=r"$\theta_1$ (mg/L)", ylabel=r"$\theta_2$ (day$^{-1}$)",
               title="Centers aligned: shape comparison" if centered else "Actual estimates and regions")
        ax.legend(fontsize=10)
    fig.tight_layout()
    return fig

fig = plot_sequential_regions(results)
helper.save_figure(fig, "bod-sequential-ellipses")

[helper] wrote media/figures/bod-sequential-ellipses.png and .pdf

A control: was the design step worth it?¶

A shrinking ellipse is not by itself evidence that the design helped. Any seventh observation adds information. The question is whether $t^\star$ adds more than a time we could have picked without optimizing.

Two controls, both inside the same 0.5--5 day window:

$t=0.5$ days, the time the dense grid identified as the worst in the window; and
$t=5$ days, the naive “measure at the end” choice.

Each control uses its own noise draw, so the three augmented data sets differ in the new point only.

grid_worst_time = grid[np.argmin(grid_logdet)]
print(f"Worst time on the dense grid: {grid_worst_time:.4f} days")

candidates = {
    f"optimal t* = {optimal_time:.3f}": optimal_time,
    f"worst on grid, t = {grid_worst_time:.3f}": grid_worst_time,
    "end of window, t = 5.0": 5.0,
}

# The prediction the design step makes, before any data are collected. It depends
# only on the sensitivity vector and the prior, not on the value that comes back.
predicted = []
for label, t_new in candidates.items():
    fim = updated_fim(t_new)
    predicted_cov = np.linalg.inv(fim)
    predicted.append(
        {
            "design": label,
            "predicted det(V)": np.linalg.det(predicted_cov),
            "predicted reduction": 1 - np.linalg.det(predicted_cov) / np.linalg.det(covariance),
        }
    )

pd.DataFrame(predicted)

Worst time on the dense grid: 0.5000 days

# Now take each measurement and re-estimate. The optimal row reuses the
# measurement already taken above, so this table and the one before it agree;
# each control draws its own noise from the same stream.
realized_covariance = {}
before = summarize("before: 6 observations", theta, covariance)
before["reduction in det(V)"] = 0.0
realized = [before]

for label, t_new in candidates.items():
    if t_new == optimal_time:
        y_new = new_observation
    else:
        y_new = run_digital_experiment(t_new, loop_rng)

    trial = pd.concat(
        [data, pd.DataFrame({"time_days": [t_new], "y": [y_new]})], ignore_index=True
    )
    _, theta_trial, covariance_trial = estimate_parameters(trial)
    realized_covariance[label] = covariance_trial

    row = summarize(label, theta_trial, covariance_trial)
    row["reduction in det(V)"] = 1 - row["det(V)"] / np.linalg.det(covariance)
    realized.append(row)

pd.DataFrame(realized)

One realization is not evidence¶

The table above compares three single draws. That is a weak comparison: $\mathbf{V}_\theta$ from the Gauss–Newton calculation depends on sensitivities at the refitted parameters as well as on where the measurement was placed, so a lucky draw at the worst time can beat an unlucky draw at $t^\star$ .

The Fisher information prediction has no such noise in it. To check that the prediction is the right summary of what actually happens, repeat each placement over many noise draws and average the realized $\det\mathbf{V}_\theta$ .

All three placements share the same standardized draws -- common random numbers -- so the comparison is between times, not between lucky and unlucky noise.

REPLICATES = 50
standardized_draws = np.random.default_rng(7).normal(0.0, 1.0, REPLICATES)

study = []
for label, t_new in candidates.items():
    determinants = []
    for z in standardized_draws:
        y_new = float(bod_response(t_new, THETA_TRUE) + MEASUREMENT_STD * z)
        trial = pd.concat(
            [data, pd.DataFrame({"time_days": [t_new], "y": [y_new]})], ignore_index=True
        )
        determinants.append(np.linalg.det(estimate_parameters(trial)[2]))

    determinants = np.array(determinants)
    study.append(
        {
            "design": label,
            "mean det(V)": determinants.mean(),
            "mean reduction": 1 - determinants.mean() / np.linalg.det(covariance),
            "predicted reduction": 1
            - np.linalg.det(np.linalg.inv(updated_fim(t_new))) / np.linalg.det(covariance),
        }
    )

pd.DataFrame(study)

def confidence_ellipse(center, cov, n_points=200):
    """Points on the 95% confidence ellipse of a 2x2 covariance matrix."""
    eigenvalues, eigenvectors = np.linalg.eigh(cov)
    angle = np.linspace(0, 2 * np.pi, n_points)
    unit_circle = np.vstack([np.cos(angle), np.sin(angle)])
    radii = np.sqrt(CHI2_95_2DOF * eigenvalues)
    return (eigenvectors @ (radii[:, None] * unit_circle)).T + np.array(center)


# All three ellipses are drawn around the six-observation estimate. Each fit has
# its own center, but overlaying them there would mix a change in shape with a
# change in location; only the shape is what the design step controls.
center = [theta["asymptote"], theta["rate_constant"]]

optimal_label, worst_label = list(candidates)[0], list(candidates)[1]
ellipses = [
    (covariance, "-", "before: 6 observations"),
    (realized_covariance[worst_label], ":", f"+ {worst_label}"),
    (realized_covariance[optimal_label], "--", f"+ {optimal_label}"),
]

for cov, style, label in ellipses:
    points = confidence_ellipse(center, cov)
    plt.plot(points[:, 0], points[:, 1], linestyle=style, label=label)

plt.scatter(*center, color="black", zorder=3)
plt.xlabel(r"$\theta_1$ (asymptote)")
plt.ylabel(r"$\theta_2$ (rate constant)")
plt.title("95% confidence regions")
plt.grid(True)
plt.legend(fontsize=9)
plt.show()

The D-optimal design minimizes the determinant of the predicted covariance, hence the area of this two-parameter ellipse. It does not specifically minimize the longest axis (the E-optimal criterion). Adding positive-semidefinite information at the same nominal parameters cannot increase variance in any fixed direction; re-estimating the parameters after collecting data can also change the sensitivities and rotate the ellipse.

Extensions and further reading¶

The BOD calculation holds the response model, noise distribution, and known noise scale fixed. The following extensions change specific parts of that formulation.

Extension	What changes	Starting point
Errors in variables	Model uncertainty in measured inputs as well as outputs; latent true inputs can become estimation variables.	Britt & Luecke (1973), DOI; Seber & Wild, §1.4
Correlated errors	Use a joint likelihood/covariance; the scalar error suffix does not encode off-diagonal terms.	White & Domowitz (1984), Nonlinear Regression with Dependent Observations; Bates & Watts, §3.8; Seber & Wild, Chapter 6
Bootstrap refinements	Choose residual, case, bias-corrected, or other intervals to match the sampling assumptions; assess simulation error.	Wehrens et al. (2000), DOI; Boos & Stefanski (2013), §11.5
Small-sample likelihood calibration	Bartlett-type corrections can reduce chi-squared calibration error under their model-specific regularity and expansion conditions.	Cordeiro & Cribari-Neto (2014), §2.4, pp. 20--23
Profile likelihood	Refit nuisance parameters to inspect asymmetry and one-sided or unbounded intervals.	Raue et al. (2009), DOI
Regularized estimation	Add a penalty and examine bias as well as variance; penalized curvature is not automatically frequentist sampling information.	Hoerl & Kennard (1970), Theorem 4.3, p. 62, DOI
Sparse estimation and priors	An $L_1$ penalty can set coefficients to zero; the choice depends on the model and sparsity assumptions.	Tibshirani (1996), §5 and §11; Park & Casella (2008), DOI
Uncertain parameters in design	Average utility over parameter uncertainty or protect against a specified uncertainty set.	Chaloner & Verdinelli (1995), article; Seber & Wild, §5.13.2
Continuous effort	Optimize nonnegative weights over candidate experiments, then convert weights into a finite allocation.	Pukelsheim (2006), §9.4 and §12.4--12.5; Kiefer & Wolfowitz (1960), DOI
Model discrimination	Choose measurements that distinguish competing response models, rather than only refining parameters of one model.	Buzzi-Ferraris & Forzatti (1983), DOI; Seber & Wild, §5.13.5
Multiresponse estimation	Estimate a common unknown error covariance for several outputs; profiling it can yield a determinant-based fitting objective.	Bates & Watts, Chapter 4, pp. 134--140
Black-box information	Approximate information from likelihood evaluations when derivative access is limited.	Spall (2005), DOI
Dynamic applications	Include trajectories, constraints, and hardware validation in model building and design.	Franceschini & Macchietto (2008), DOI; Wilson et al. (2014); Bandara et al. (2009)

For a Gaussian prior with precision $\mathbf P$ and mean $\boldsymbol\mu$ , MAP estimation adds $\tfrac12(\boldsymbol\theta-\boldsymbol\mu)^T\mathbf P(\boldsymbol\theta-\boldsymbol\mu)$ to the negative log likelihood. Its Hessian contributes $\mathbf P$ to local posterior curvature. Ridge theory can establish a bias--variance tradeoff in specified linear models; it does not guarantee a large variance reduction for an arbitrary nonlinear BOD fit. An $L_1$ penalty corresponds to independent Laplace priors for MAP estimation, while Bayesian posterior sampling is a different operation.

The ParmEst tutorial and Pyomo.DoE tutorial retain the larger dynamic-model implementations. Classical factorial designs, blocking, and response surfaces address additional experimental-design questions; see Dean, Voss & Draguljić (2017), Chapters on factorial, block, and response-surface designs. An information-based design for one fixed model does not replace those methods.

Carry remaining uncertainty into decisions¶

Stopping the measurement loop depends on prediction needs and experimental budget. Bootstrap parameter draws can support scenario construction if their sampling interpretation is suitable for the decision task. A confidence region alone supplies no scenario probabilities. Stochastic programming develops decisions under specified scenarios and probabilities.

Key takeaways¶

Parameter estimation turns measurements into fitted parameters and an uncertainty description.
The inverse Gauss–Newton covariance supplies local prior information for the next design.
Experimental design asks where a new measurement adds information not already supplied by the prior.
A design is only a prediction. Closing the loop -- taking the measurement and re-estimating -- is what tests it, and a digital experiment makes that possible when the real apparatus is gone.
Check the design against a control. That a seventh observation shrinks the confidence region proves nothing on its own; the comparison that matters is against a measurement placed without optimizing.
A single realization is a weak comparison, because the estimated covariance changes with the refitted parameters. Average over draws, or compare the noise-free Fisher information prediction.
FIM designs are local: different nominal parameters, noise models, bounds, or priors can change the answer.
Small examples permit independent checks. Use them before trusting a larger dynamic design.

References¶

Wehrens, R., Putter, H., & Buydens, L. M. C. (2000). The bootstrap: a tutorial. Chemometrics and Intelligent Laboratory Systems, 54(1), 35–52. 10.1016/s0169-7439(00)00102-7
Laky, D. J., Lilonfe, S., Martin, S. B., Klise, K. A., Nicholson, B. L., Siirola, J. D., & Dowling, A. W. (2026). Optimal experimental design using eigenvalue-based criteria with Pyomo.DoE. Digital Discovery, 5(9), 3780–3803. 10.1039/d6dd00160b
Britt, H. I., & Luecke, R. H. (1973). The Estimation of Parameters in Nonlinear, Implicit Models. Technometrics, 15(2), 233–247. 10.1080/00401706.1973.10489037
Raue, A., Kreutz, C., Maiwald, T., Bachmann, J., Schilling, M., Klingmüller, U., & Timmer, J. (2009). Structural and practical identifiability analysis of partially observed dynamical models by exploiting the profile likelihood. Bioinformatics, 25(15), 1923–1929. 10.1093/bioinformatics/btp358
Hoerl, A. E., & Kennard, R. W. (1970). Ridge Regression: Biased Estimation for Nonorthogonal Problems. Technometrics, 12(1), 55–67. 10.1080/00401706.1970.10488634
Kiefer, J., & Wolfowitz, J. (1960). The Equivalence of Two Extremum Problems. Canadian Journal of Mathematics, 12, 363–366. 10.4153/cjm-1960-030-4
Buzzi-Ferraris, G., & Forzatti, P. (1983). A new sequential experimental design procedure for discriminating among rival models. Chemical Engineering Science, 38(2), 225–232. 10.1016/0009-2509(83)85004-0
Spall, J. C. (2005). Monte Carlo Computation of the Fisher Information Matrix in Nonstandard Settings. Journal of Computational and Graphical Statistics, 14(4), 889–909. 10.1198/106186005x78800
Franceschini, G., & Macchietto, S. (2008). Model-based design of experiments for parameter precision: State of the art. Chemical Engineering Science, 63(19), 4846–4872. 10.1016/j.ces.2007.11.034