mixle.doe.calibrate module

Kennedy-O’Hagan calibration: infer a simulator’s parameters with an explicit model-discrepancy term.

Field data rarely equals the simulator even at the true parameters – there is model-form error. Fitting parameters by plain least squares absorbs that bias and gives wrong (over-confident) parameters. The Kennedy-O’Hagan model writes y(x) = eta(x, theta) + delta(x) + noise with delta a GP discrepancy, and infers theta and delta jointly, so the parameters are not contaminated by the bias.

calibrate(simulator, x, y, theta0, *, discrepancy=True, discrepancy_lengthscale=None, seed=0, max_iter=1000)[source]

Calibrate simulator(x, theta) to field data (x, y) with a GP discrepancy term.

Maximizes the marginal likelihood of the residual r(theta) = y - eta(x, theta) under a GP + noise model, over theta and the discrepancy amplitude + noise. discrepancy=False drops the GP (plain nonlinear least squares) – useful to show the bias the discrepancy removes.

The discrepancy correlation length is fixed (discrepancy_lengthscale, default 10% of the input domain) rather than fitted: this is the standard resolution of the Kennedy-O’Hagan theta/delta identifiability problem – a short discrepancy length forces the GP to model only local model-form error, leaving the smooth global trend to the parametric simulator so theta stays identifiable. Set it to the scale of model error you expect.

Parameters:

  • simulator (Callable[[ndarray, ndarray], ndarray]) – eta(x, theta) -> predictions (vectorized over the rows of x).

  • x (ndarray) – field inputs and observations. x and y must have the same number of rows, and x, y, theta0, and simulator(x, theta0) must all be finite.

  • y (ndarray) – field inputs and observations. x and y must have the same number of rows, and x, y, theta0, and simulator(x, theta0) must all be finite.

  • theta0 (Sequence[float]) – initial calibration parameters (its length sets the parameter count). There must be more observations than free theta parameters, and the fitted simulator sensitivity matrix must have full parameter-column rank.

  • discrepancy (bool) – include the GP discrepancy term (the Kennedy-O’Hagan model).

  • discrepancy_lengthscale (float | None) – fixed GP correlation length; None (default) uses 10% of the input domain. If given explicitly it must be positive – unlike None, 0 is not silently treated as “use the default”: a caller-supplied non-positive value is a contract violation, not a request for the default.

  • seed (int) – seeds the small, reproducible random perturbation used to initialize the optimizer, so the same seed value reproduces the same fit and a different seed can escape an exact tie or a degenerate flat ridge at the unperturbed starting point.

  • max_iter (int) – positive optimizer iteration budget. Reused, unmodified, for the internal no-discrepancy warm-start fit when discrepancy=True.

Raises:

ValueError – on an x/y row-count mismatch, non-finite x/y/theta0/ simulator(x, theta0), a simulator(x, theta0) shape that does not match y, a non-positive discrepancy_lengthscale or max_iter, fewer observations than theta parameters, rank-deficient simulator sensitivity, or if the optimizer fails to converge to a finite objective.

Return type:

KOCalibration

exception CalibrationIdentifiabilityError(*, rank, n_parameters, singular_values, non_identifiable_directions)[source]

Bases: ValueError

Raised when simulator sensitivity cannot identify every calibration-parameter direction.

Parameters:

  • rank (int)

  • n_parameters (int)

  • singular_values (np.ndarray)

  • non_identifiable_directions (np.ndarray)

Return type:

None

class KOCalibration(theta, ls, amp, noise, simulator, x, y, theta_standard_error=None, *, sensitivity_rank=None, sensitivity_singular_values=None, sensitivity_condition_number=None)[source]

Bases: object

Result of calibrate(): a POINT MLE of the fitted parameters (not a posterior), the discrepancy GP, and a calibrated predictor.

theta is the maximizer of the (marginal) likelihood, not a full Bayesian inference of it – there is no prior and no posterior distribution over theta here. theta_standard_error is an asymptotic (Laplace / observed-Fisher-information) standard error for each theta component, from the inverse Hessian of the negative log-likelihood at the MLE; theta_ci_low/ci_high are the corresponding two-sided 95% Wald intervals (theta +- 1.96 * theta_standard_error). This is a local, asymptotic approximation (it can be optimistic for small n, strong parameter correlation, or a non-quadratic likelihood near the optimum) – not a substitute for a real posterior when that matters. All three are arrays of nan when the Hessian at the optimum is not numerically invertible (or not positive-definite), which the point estimate itself does not depend on. sensitivity_rank, sensitivity_singular_values, and sensitivity_condition_number certify local structural identifiability; a rank-deficient fit raises CalibrationIdentifiabilityError instead of returning an arbitrary point.

All fitted arrays are privately copied and frozen. Array properties return detached copies, so neither caller-owned inputs nor a returned view can retroactively change the fitted residual, interval, or predictor.

Parameters:

  • sensitivity_rank (int | None)

  • sensitivity_singular_values (np.ndarray | None)

  • sensitivity_condition_number (float | None)

property theta: ndarray

Detached fitted parameter vector.

predict(x_new, *, with_discrepancy=True)[source]

Calibrated prediction at x_new: simulator at the fitted theta, plus the GP discrepancy (the bias-corrected estimate of reality) unless with_discrepancy=False (the pure simulator).

Parameters:
Return type:

ndarray