Model calibration: general concepts

ACBICI is a suite of Python routines that can be used to calibrate (complex) models, finding the optimal probability distributions of their unknown parameters, given some experimental data and the model itself. The library design favors clarity and simplicity, and thus users can adapt it to fit their own needs with minor modifications.

Before describing how to set up a calibration process, let us introduce a few concepts.

Models

A model is a mathematical expression, equation, or code that embodies the behavior of a system, a person, an object, etc. More formally, a model defines a mapping between an input space and a target space, usually with the goal of predicting or analyzing the behavior of the modeled system.

The input and target spaces can be, in general, abstract sets. In ACBICI, however, we restrict the input space to \(\mathbb{R}^d\), \(d\ge 1\), and the target space to \(\mathbb{R}\). That is, at this moment, the library accepts multidimensional inputs in a real vector space and only scalar outputs.

Often, models depend on parameters which we refer to as \(\theta\). These parameters belong to yet another space which in ACBICI can only be \(\mathbb{R}^p\) with \(p\ge1\). Both parameters and input variables determine the value of a model output but their roles are different. Input variables have values that are known, at least to a certain extent, coming from measurements, controls, pools, actuators, etc. In contrast, parameters are artificial variables introduced in the model to make it predictive, but need not be measurable. In addition, once the parameters are set in a model, they remain fixed for all inputs and predictions. We say that the model is calibrated and the accuracy of its predictions depends both on the goodness of the model itself as well as in the quality of the selected parameters.

The point of view of Bayesian calibration – and thus, ACBICI – pivots on the assumption that parameters are always random variables. Their probability distribution is often unknown and the goal of inference is to determine its optimal form, given the experimental evidence available and possibly some prior knowledge. In Bayesian statistics, the criterion for optimality comes from Bayes rule, which combines in a precise way the information from the experimental data, the prior information, and the model.

Calibrating a model

As indicated before, the goal of ACBICI is to calibrate models. To be specific, the model is denoted as \(g\), and it can refer to a simple mathematical function that can be written in analytical form or a complex code whose inner workings are not fully known. In any case, a model is a mapping \(\mathbb{R}^d\times \mathbb{R}^p \ni (x,\theta)\mapsto y=g(x;\theta)\in \mathbb{R}\), where \(x\) is referred to as the input variables and \(\theta\) are the parameters.

The goal of calibration is to find the optimal value of the parameters of a given model \(g\), when there is some experimental data in the form of pairs \((x_i,y_i)\in \mathbb{R}^d\times \mathbb{R}, i=1,2,...N\). Different notions of optimality produce different calibration results. In particular, in ACBICI, the library looks for the probability distribution for \(\theta\) that best combines the model, the prior information we have about this distribution, and the experimental data.

To calibrate a model in ACBICI three things are needed:

1. A model, that is, a Python function that can be called with some values for \(x\) and \(\theta\) and can calculate \(y=g(x;\theta)\). In some cases it is enough to have the results of evaluating the model for a rich enough set of pairs \((x_i,\theta_i)\). More about this later.

2. Experimental data stored in a data file containing pairs \((x_i,y_i)\).

3. Prior information about the probability distribution of the parameters \(\theta\). This information can be very precise or vague, but due to the Bayesian structure of the code some prior information is always required. In ACBICI, it is enough to select the type of prior information that one would like to consider since the library will provide helper functions to implement them.

To calibrate a model with ACBICI, thus, the user must prepare the three things enumerated above. In addition, the user must make three decisions that will determine the type of calibration that will be performed:

1. All experimental data has some error, irrespective of its origin. In ACBICI, experimental data is assumed to have a Gaussian probability distribution with zero mean. The user must decide whether the covariance of this error is known, and thus provide it, or it is unknown and should be estimated.

2. The user must assess whether the model we would like to calibrate is considered to be costly to run or now.

3. Finally, since all models are incomplete to a certain extent, the user might decide to accept the proposed model as good enough to represent the data or would like to estimate, in addition to all calibration tasks, the discrepancy of the model with the experimental data.

Let us provide some additional explanations regarding the two last topics.

Expensive or inexpensive models

To understand the need for this choice, we recall that in all kinds of Bayesian inference the goal is to obtain a posterior probability distribution which, in most cases, has no analytical expression. Rather, instead of a closed-form formula, in Bayesian calibration we obtain a (large) sample of points whose probability distribution is precisely the one we are looking for. To obtain it, a Markov chain Monte Carlo (MCMC) method needs to be employed. This type of methods requires evaluating the model hundreds or even thousands times. If the model that needs to be calibrated takes a long time to be evaluated, it might be prohibitively expensive to run MCMC with it and we then consider it as expensive.

In ACBICI, an expensive model is replaced by a Gaussian process that acts as surrogate. This meta-model is much faster to evaluate than the original model but its hyper-parameters need to be calibrated alongside the original parameters \(\theta\). Everything comes with a price: either you use the original model in the MCMC computations or you evaluate it a few times to calibrate the surrogate. For very expensive models it might be mandatory to use a meta-model at the expense of losing some accuracy in the overall calibration process. In these situations, the user has to supply a file with data triplets \((x_i,\theta_i,g(x_i;\theta_i))\) that are obtained evaluating the model \(g\) at carefully selected pairs \((x_i,\theta_i)\). These pairs must cover the input and parameter space well enough. More details will be provided later.

Discrepancy error

Every model, irrespective of its sophistication, ignores deliberatively or not details of the true system it represents, be it for simplicity or lack of knowledge. As a result, all models are inexact or, in other words, “All models are wrong, but some are useful” (G.E.P Box).

When performing the calibration of a model, we can ignore the limitations of the latter and insist in obtaining the best possible distribution for each of its parameters. It is also possible, however, to model also the discrepancy between the experimental results and the model predictions. This discrepancy error in ACBICI is also represented by a Gaussian process that, if present, provides additional hyperparameters that need to be calibrated, thus increasing the cost of the process.

In summary, depending on the calibration choices we might need to build a surrogate and the discrepancy model. The following table summarizes all the calibration types in ACBICI:

Calibration types

Name	Expensive model	Discrepancy surrogate	Hyperparameters
Type A	NO	NO	None
Type B	YES	NO	3 (meta-model)
Type C	NO	YES	2 (discrepancy)
Type D	YES	YES	3 (metamodel) + 2 (discrepancy)

Let us note that the complexity of calibrating a model grows with the number of parameters and hyperparameters. In fact, the larger this number is, the more costly is to evaluate the calibration model and the larger the number of MCMC steps to be performed. This comment should caution the user from using calibration strategies C and D. On the other hand, the knowledge of the discrepancy error can be a useful quantity to ascertain, in a global way, the adequacy of a model to reproduce experimental data.

Identifiability analysis

Experimental data may not contain enough information to distinguish every model parameter, hyperparameter, or inferred experimental-error term. In that case, the posterior can remain close to its prior, and point estimates or credible intervals should not be interpreted as evidence that the parameter is well determined. ACBICI therefore provides an identifiability analysis that is run after calibrate() and uses the posterior samples saved by the calibration.

The first diagnostic is the prior-to-posterior compression factor. For every calibrated quantity, ACBICI computes

\[c = 1 - \frac{\operatorname{std}(\mathrm{posterior})} {\operatorname{std}(\mathrm{prior})}.\]

A value near one indicates substantial contraction: the posterior is much narrower than the prior and the data have strongly constrained that quantity. A value near zero indicates little contraction, so the data have added little information beyond the prior. A negative value is possible when the posterior is wider than the prior and should prompt further investigation. This diagnostic requires a finite, nonzero prior standard deviation; otherwise its compression factor is reported as unavailable.

The analysis appends the compression factors to the calibration log and writes identifiability.png. The figure overlays the prior density with the posterior histogram and density estimate for each model parameter, GP hyperparameter, and inferred experimental-error parameter, where applicable. It is useful to read the numerical factor and the overlaid distributions together: a narrow posterior concentrated near a prior boundary, for example, may need additional modelling scrutiny even when its compression factor is large.

The second diagnostic assesses sensitivity to the chosen priors. By default, ACBICI rescales each prior spread by factors of 0.9 and 1.1 and uses importance reweighting of the existing posterior samples to estimate the corresponding changes in posterior means and standard deviations. It records these changes, along with the effective sample size (ESS) of the reweighting, in the log and writes sensitivity.png. This is fast because it does not require new MCMC runs. However, results for a scaling factor whose reweighting ESS is below 10 percent should be treated as unreliable. The standalone analysis functions can instead run fresh MCMC chains for the scaled priors when a more robust, but more expensive, sensitivity calculation is required.

GUI

ACBICI includes a graphical user interface (GUI) in the gui directory. It is a local browser-based workbench for defining, running, and reviewing a calibration without having to write the driver script by hand. It is intended to make the normal ACBICI workflow accessible while keeping the resulting analysis reproducible: every analysis configured in the GUI is converted into an ordinary Python script that can be inspected, edited, and run independently.

Start the GUI from the repository root with:

python gui/acbici_gui.py

By default, the command opens a browser at http://127.0.0.1:8766. The server binds only to the local machine, so it does not expose the project files on the network. Use --port to choose a different port, or --no-browser when working on a machine where the browser must be opened manually. Stop the GUI server with Ctrl+C in the terminal from which it was launched.

The interface is organized into Setup, Run, Results, Report, and Manual views. The Setup view is where an analysis is defined. First, choose one of the four calibration strategies described above: classical Type A, expensive-model Type B, discrepancy Type C, or Kennedy–O’Hagan Type D. For the GP-based strategies, the view also exposes the covariance-kernel selection. The current choice is displayed beside the run name, making it easy to confirm the kind of analysis before it is started.

The model panel lets the user browse to a Python file containing an ACBICImodel subclass and load it into the workbench. Its source is displayed in an editable text area, and the Save edits control writes any changes back to the selected model file. After loading, the GUI reads the model definition and presents one editor per model parameter. Each editor can select any ACBICI prior distribution and its numerical arguments, allowing the priors declared in the model to be retained or overridden. This is particularly useful for trying prior choices while preserving a model implementation in a separate file.

Experimental observations are selected from a data file using the local file browser. The experimental-error setting can retain the library default, use a known standard deviation, or infer an unknown error using a chosen prior. For Type B and Type D analyses, a synthetic-data panel also accepts precomputed surrogate training data. If it is left empty, the user can set the number of synthetic points to be generated by Latin hypercube sampling. GP hyperparameter-prior editors are shown for calibrator types that use them; the library defaults remain in effect unless an individual override is enabled.

The calibration panel selects between ensemble MCMC (emcee), VBMC (pyvbmc), and the internal diagonal-covariance VBMC implementation. MCMC controls include the number of steps, burn-in fraction, thinning interval, and number of walkers. The plots-and-diagnostics panel selects which postprocessing figures should be produced and can request the identifiability analysis, including its default reweighting-based sensitivity diagnostic. These options correspond directly to arguments in the generated calibration script.

Select Generate script after completing the setup. The GUI writes a self-contained script to <run name>.out/run_<run name>.py and shows its contents in the generated-script panel. The script imports the selected model, applies the selected priors and error treatment, loads the data, creates the appropriate calibrator, runs the chosen inference method, and produces the requested plots. Because paths are recorded explicitly, the script is a useful record of the analysis and can be rerun outside the GUI. The Load previous script and Recent analyses controls reconstruct the GUI settings from scripts previously generated by the workbench, so an earlier run can be used as a starting point for a modified analysis.

The Run view launches the generated script as a background process. It reports the current status, elapsed time, command, exit code, and live solver output. MCMC progress is displayed when it is available, and the running job can be stopped from this view. Once the process has completed, Results collects the output figures in one place, including applicable model-versus-data, discrepancy, prior/posterior, trace, autocorrelation, split-\(\hat{R}\), corner, correlation, prediction-error, identifiability, and sensitivity plots. Each figure is labelled with a short explanation of its diagnostic purpose. Finally, Report displays the acbici.log file produced by the run, while the Manual view provides the built documentation without leaving the workbench.

Outputs from the calibration

Once a model has been calibrated, there is a wealth of information that can be obtained from the posterior probability distribution. In fact, ACBICI computes the joint probability distribution of all the parameters and hyperparameters so any statistic can be, in principle, obtained from the former.

• Parameter statistics, including the mean, variance, and maximum a posteriori estimate (MAP).

• Corner plots: the outcome of the calibration is a multi-dimensional probability distribution for all the parameters of the model and hyperparameters of the meta-model. To better understand the correlations between these variables, a corner plot is depicted that includes all the one- and two-dimensional projections of the sample.

• Correlation plots: when calibrating a model with multiple parameters and/or hyperparameters, it might be interesting to assess the strength of the correlation among variables. A correlation plot shows with a color scheme the value of the Pearson correlation parameter for each pair of calibrated variables. Highly correlated variables might indicate that they are not independent and thus might lead to issues in the calibration.

• Comparison plots: a plot is generated that depicts the experimental data and the output of the calibrated model, as a curve. This plot is not shown when the dimension of the input space is larger than one.

• Pearson plot: The correlation between all parameters and hyperparameters is shown as a heat map. Large correlations (+1 or -1) in off-diagonal blocks are indication that the corresponding parameter are not independent, signaling potential errors.

• Prior/posterior plots: This plot shows the prior distribution functions for each parameter and hyperparameter next to its posterior probability densities. One should pay attention to situations in which the posterior is very similar to the prior in which case the sampled data might not have provided enough information to update the former.

• Prediction plot: This plot shows a histogram of the errors made by the calibrated model as compared with the experimental data. Also the mean error in the predictions is indicated.

• Trace plot: This plot shows, for each parameter and hyperparameter, its value along the MCMC chain. Users should pay attention to results that are not well-mixed, meaning that the variables do not explore their potential values. In these situations, the MCMC algorithm may have failed to identify regions of large likelihood and thus failing to sample the posterior.

Convergence

MCMC diagnostics

The calibration report includes several standard Markov chain Monte Carlo (MCMC) diagnostics that help assess whether the posterior samples are reliable.

• Effective Sample Size (ESS): The report provides the ESS for each parameter as well as a mean ESS value. ESS estimates how many effectively independent samples the chain contains after accounting for autocorrelation. High ESS values generally indicate that the posterior summaries (means, credible intervals, etc.) are well estimated.

• split-:math:`hat{R}` convergence diagnostic: For each parameter we report the split-\(\hat{R}\) statistic. Values very close to \(1\) indicate that multiple chains have mixed well and are consistent with sampling from the same stationary distribution.

• Autocorrelation time (:math:`tau`): The mean integrated autocorrelation time \(\tau\) describes how quickly the sampler explores the posterior distribution. Shorter \(\tau\) implies better sampling efficiency and efficient exploration. We also use \(\tau\) as a practical convergence / run-length diagnostic: the chain should be long enough compared to \(\tau\) to ensure stable estimates.

  More details are available in the emcee autocorrelation tutorial and in Goodman & Weare (2010), Ensemble samplers with affine invariance.

Note

In addition to the scalar diagnostics above, we generate autocorrelation plots and trace plots.

• Autocorrelation plots: should show that autocorrelation decays with lag and that the estimated integrated autocorrelation time stabilizes as the chain length increases.

• Trace plots: show each parameter value versus step. Trace plots are used as a quick visual diagnostic to assess mixing, stationarity after burn-in, and the presence of long-term trends or poor exploration.

• Gelman-Rubin statistics plots: summarize the split-\(\hat{R}\) convergence diagnostics for all calibrated and surrogate-related parameters. Values close to one and below the indicated threshold lines, indicate that no parameter exhibits problematic convergence behavior and that the MCMC chains have mixed satisfactorily.

Warning

If you see an error like:

emcee.autocorr.AutocorrError: The chain is shorter than 50 times the integrated autocorrelation time for 1 parameter(s).
Use this estimate with caution and run a longer chain! N / 50 = 22;
tau: [22.14675566 23.67994754]

the chain is too short relative to the estimated \(\tau\). In that case, rerun the calibration with at least:

n_mcmc = max(tau) * len(p) * (1 + burnFraction) * nwalkers

For example:

24 * 2 * 1.2 * 32 = 1843

(rounded to full digits) gives \(n_\mathrm{mcmc}=1843\) steps. Repeat until the run is long enough that the autocorrelation-time estimate is valid and stable.

VBMC diagnostics

For variational inference with VBMC, convergence is assessed using diagnostics that are specific to optimization-based posterior approximation, rather than MCMC chain behavior. VBMC iteratively refines a variational posterior and monitors several quantities to determine when the approximation has stabilized.

• Evidence Lower Bound (ELBO): During optimization, VBMC tracks the evidence lower bound (ELBO), which serves as the objective function for variational inference. Convergence is indicated when the ELBO reaches a plateau and subsequent changes become small relative to the estimated ELBO uncertainty (standard deviation). A stabilized ELBO suggests that further optimization is unlikely to substantially improve the posterior approximation.

• Symmetrized KL divergence (sKL): VBMC monitors the symmetrized Kullback–Leibler (sKL) divergence between variational posteriors from successive iterations. As optimization progresses, the sKL should decrease toward zero, indicating that the variational posterior is no longer changing appreciably and has reached a stable solution.

• Reliability (convergence) index (:math:`r_mathrm{index}`): The ELBO stability and inter-iteration sKL are combined into a single reliability (convergence) index, \(r_\mathrm{index}\). In :texttt{PyVBMC}, inference is terminated only once \(r_\mathrm{index} < 1\) and remains below this threshold for a prescribed stability window. This criterion ensures that both the optimization objective and the variational posterior have converged.

Note

VBMC diagnostics play a role analogous to ESS, \(\hat{R}\), and autocorrelation time in MCMC. While MCMC diagnostics assess sampling efficiency and mixing, VBMC diagnostics assess the stability and accuracy of the variational approximation during optimization.

In practice, stabilized ELBO values, near-zero inter-iteration sKL, and a low \(r_\mathrm{index}\) collectively indicate robust convergence for VBMC. In well-behaved analyses, the resulting marginal posterior distributions should exhibit clear peaks for most parameters, signaling that the algorithm has accurately approximated the dominant modes of the posterior distribution.

Multi-Output Calibration

When not explicitly stated, ACBICI operates on single-output calibration problems. Nevertheless, the framework naturally extends to multi-output calibration.

With the exception of the comparison plot, all features described above are also available in the multi-output setting.

Multi-output calibration in ACBICI is performed by augmenting the input space with a discrete output (task) index, which reduces the multi-output problem to an equivalent single-output formulation. The augmented input space is

\[(x,i) = r \in \mathbb{R}^d \times \{0,\ldots,n_{\text{tasks}}-1\},\]

where \(n_{\text{tasks}}\) denotes the number of calibration outputs.

Each output is associated with a task-specific model \(m_i : \mathbb{R}^d \times \mathbb{R}^p \to \mathbb{R}\), for \(i = 0,\ldots,n_{\text{tasks}}-1\). These are combined into a single effective model on the augmented space,

\[m : \mathbb{R}^d \times \{0,\ldots,n_{\text{tasks}}-1\} \times \mathbb{R}^p \to \mathbb{R},\]

defined as

\[m(x,i,\theta) = \sum_{j=0}^{n_{\text{tasks}}-1} m_j(x,\theta)\,\delta_{ij},\]

where \(\delta_{ij}\) is the Kronecker delta. For a fixed task index \(i\), the corresponding restriction of \(m\) coincides with the original model \(m_i\).

In practice, vector-valued experimental and synthetic datasets are mapped to the augmented input space by replicating each input once per task and assigning the corresponding task index. The multi-output observations are then stacked into a single scalar response vector, enabling calibration with a joint likelihood.

When surrogate models are used (types B, C, and D), this augmentation requires a task-dependent covariance function. In ACBICI, a similarity kernel is employed, where a Matérn 3/2 kernel is used when the task indices match (\(i=j\)), and the covariance is zero otherwise:

\[\begin{split}k((x, i), (x', j), \lambda, \beta) = \begin{cases} k_{\text{m32}}((x,x'), \lambda, \beta), & \text{if } i = j, \\ 0, & \text{if } i \neq j. \end{cases}\end{split}\]

Here, \(k_{\text{m32}}\) denotes the Matérn 3/2 kernel parameterized by \(\lambda\) and \(\beta\).

In multi-output calibration, nondimensionalization and output scaling are especially important. Since all outputs are weighted equally in this formulation, they must be normalized to comparable scales to ensure balanced and meaningful calibration results.