Optimal experimental design: turn profile_likelihood's "practically unidentifiable" verdict into the experiment that fixes it, reusing the assemble_fisher_hessian already built for gntr

[wshlavacek]

job_type = profile_likelihood (#446) ends by telling a user that a parameter is practically unidentifiable — a finite lower crossing, an open CI at the bound, a flat profile. That is a diagnosis with no prescription. The question an experimentalist immediately asks is "then what experiment should I run?", and PyBNF has no answer.

Optimal experimental design (OED) answers it, and — this is the reason to open the issue now — the expensive machinery already exists.

The Fisher information is already assembled

assemble_fisher_hessian is a standalone public API (pybnf/gradient/__init__.py), built for the EFIM trust-region optimizer gntr (#481). It assembles the expected-Fisher / Gauss-Newton information in the same point walk as the scalar gradient, and it is not a fixed-σ Gaussian special case: per gradient/__init__.py, the surface covers an estimated noise scale (#451), a log/lognormal scale (#452), per-observable trajectory transforms and normalization (#453), the asymmetric Laplace / Student-t families with mean centering (#454), and constraint penalties (#456), with per-family kappa_i location-Fisher factors.

That matrix — F(θ, ξ) for a design ξ — is the object OED optimizes. The hard, general, noise-family-aware part is done. What is missing is the loop around it.

Proposal

A design ξ is a choice over the levers PyBNF already models: which observable to measure, at which timepoints, under which condition/perturbation (the existing condition/experiment: surface), and how to allocate replicates. Score a candidate design by a scalar functional of the information matrix:

criterion	objective	reads as
D-optimal	maximize log det F	shrink the joint confidence ellipsoid's volume
A-optimal	minimize tr(F⁻¹)	shrink average parameter variance
E-optimal	maximize λ_min(F)	attack the worst direction — the sloppiest one
c-optimal	minimize cᵀF⁻¹c	target one parameter (or contrast) the PL run flagged

The natural PyBNF surface is c- and E-optimal, because they connect directly to a profile_likelihood verdict: "k_deg came back practically unidentifiable" becomes "maximize information along k_deg," which is a well-posed, single-answer question.

Two regimes, in order:

1. Local OED at θ̂ — evaluate F at the best fit. Cheap, and immediately useful. Its known weakness is that the design is only as good as θ̂; a design optimized at a bad point estimate can be worthless.
2. Robust / Bayesian OED — average the criterion over an ensemble. PyBNF already produces two suitable ensembles for free: posterior draws from am/dream/p_dream, and the _ProfileTrack points a profile likelihood run already computed and wrote to Results/profile_<name>.txt. Reusing the profile points is the cheap version and is, as far as I know, not something other tools do.

Design questions that need deciding

• Surface. A new job_type = design (a run that consumes a fit and emits a ranked design table), versus a post-fit report emitted alongside profile_likelihood_summary.txt. I lean toward the latter first — it makes OED a finding of the identifiability run rather than a separate job the user must know to ask for.
• The combinatorics are real. "Pick the best 5 of 200 candidate timepoints" is a subset-selection problem, not a continuous one. Standard resolutions: a greedy exchange algorithm (Fedorov), or the convex relaxation to continuous design weights (Kiefer) followed by rounding. Greedy is submodular-ish for D-optimality and is the pragmatic starting point.
• Enumerating the candidate space. Timepoints within the simulated window are easy. New conditions are harder — proposing a perturbation means proposing a model modification, which PyBNF has no vocabulary for. Suggested first cut: timepoints and observable selection only, within existing conditions. That is the 80% case and it sidesteps the open-ended part.
• Singular F. A structurally unidentifiable parameter makes F singular and tr(F⁻¹) infinite — which is correct information, not a numerical failure. It should be reported as "no design over this space can identify this parameter" and routed to the structural-identifiability check rather than crashing. (See the sibling issue on structural identifiability.)
• Does it need to re-simulate? F at a candidate timepoint needs the sensitivity there. If candidate points lie inside the already-simulated window, the existing sensitivity tensor (A — Gradient plumbing: capability gate + preserve the output-sensitivity tensor through net execution #447) can be read at those rows and no new solve is needed — a large cost saving worth designing for deliberately.

Why this is worth doing

It changes what the tool is for. Every other item on PyBNF's roadmap makes fitting better; this one makes the next experiment better, which is the thing that actually determines whether the model can be identified at all. It is also the natural sequel to the feature that just shipped — profile_likelihood produces exactly the input OED consumes — and it reuses assemble_fisher_hessian rather than requiring new sensitivity math.

Evaluation before committing

1. On a corpus slug with a known unidentifiable parameter, confirm the local-OED recommendation is sane — that adding the recommended timepoints measurably narrows that parameter's profile CI on a refit against synthetic data generated at θ̂.
2. Quantify the gap between the local-OED design and the robust design over the profile ensemble. If local OED is unstable across the ensemble, ship robust-only and say so.
3. Confirm assemble_fisher_hessian can be evaluated at arbitrary candidate rows without a re-solve.

Non-goal. Closed-loop / adaptive design that automatically runs the experiment. This emits a ranked recommendation; a human runs the assay.