Bayesian Calibration

Good question. Roughly speaking, Bayesian calibration estimates parameters as a "probability distribution" rather than a "single value". Ordinary least squares or optimization gives a point estimate like "This is the value that best fits the experimental data!", but Bayes returns an answer as a distribution saying "The parameter is within this range with this probability."

For example, let's say we estimate the Young's modulus of an FEM model for automotive crash analysis. Optimization would only give a single value like $E = 210\,\text{GPa}$. But Bayes would output "a distribution of $E$ with mean $210\,\text{GPa}$ and standard deviation $5\,\text{GPa}$". This uncertainty information can be directly passed to downstream reliability analysis.

More importantly, it can formally incorporate engineers' experiential knowledge (prior distribution). It narrows down parameters by utilizing both prior knowledge like "Young's modulus for this material is probably around $200 \sim 220\,\text{GPa}$" and experimental data. This is the greatest strength of Bayes.

The mathematical backbone of Bayesian calibration is Bayes' theorem. For a parameter vector $\boldsymbol{\theta}$ and observed data $\mathbf{D}$:

Let's consider the simplest case. Assuming the difference between $N$ experimental data points $d_i$ and model predictions $M(\mathbf{x}_i, \boldsymbol{\theta})$ follows independent and identically distributed normal errors:

Here, $\sigma$ is the standard deviation of the observation noise. It can be fixed as known or estimated together with $\boldsymbol{\theta}$ using Bayes. In practice, it's more robust to include $\sigma$ as an estimation target because the exact magnitude of error is often unknown.

Sharp observation! That's exactly what the Kennedy-O'Hagan framework systematically addresses.

The framework proposed by Kennedy and O'Hagan in 2001 has become the de facto standard for Bayesian calibration in CAE. The core idea is to decompose the experimental value as follows:

If you don't include $\delta$, all structural model error gets forced onto the parameters $\boldsymbol{\theta}$. For example, if the discrepancy is actually due to inadequate boundary conditions, the Young's modulus might be pushed to an extreme value to make things fit. This is called the parameter compensation problem, leading to parameter estimates that are physically meaningless.

By explicitly including $\delta(\mathbf{x})$, the "part that the model cannot reproduce" is separated, keeping parameter estimation within a physically plausible range.

$\delta(\mathbf{x})$ is learned from data as a Gaussian Process (GP). That is, the shape of $\delta$ is not predetermined but is automatically estimated via Bayesian inference from the residual pattern between experimental data and model predictions. The hyperparameters of the GP's kernel function (e.g., squared exponential kernel) are also estimated together.

Setting the prior is indeed the most debated point in Bayes. In practice, the following distinctions are often used:

• Informative Prior: Set based on material database values, past experimental results, handbook values. Example: For structural steel Young's modulus, $E \sim \mathcal{N}(210, 5^2)\,[\text{GPa}]$.
• Weakly Informative Prior: Loosely constrains the physically plausible range. Example: Friction coefficient $\mu \sim \text{Uniform}(0.05, 0.8)$.
• Non-informative Prior: Such as Jeffreys prior, used when you want to let the data speak. However, convergence of the posterior can be slow in high dimensions.

In the field, weakly informative priors are common. It's a balance of "excluding physically impossible regions while giving sufficient influence to the data."

When data is scarce, the influence of the prior is large. Conversely, with sufficient data, the posterior distribution converges to almost the same place regardless of the prior. This is called asymptotic consistency in Bayes. Therefore, in practice, you should always perform prior sensitivity analysis (checking how much results change when varying the prior).

Exactly. The normalization constant $p(\mathbf{D})$ of the posterior is a high-dimensional integral that cannot be solved analytically. That's where Markov Chain Monte Carlo (MCMC) methods come in. They generate samples directly from the posterior distribution; increasing the number of samples allows approximating the posterior to any desired accuracy.

The flow of the most basic Metropolis-Hastings algorithm is as follows:

1. Set initial value $\boldsymbol{\theta}^{(0)}$
2. Generate a candidate point $\boldsymbol{\theta}^*$ from a proposal distribution $q(\boldsymbol{\theta}^*|\boldsymbol{\theta}^{(t)})$
3. Calculate the acceptance probability:

Accept the candidate $\boldsymbol{\theta}^*$ with probability $\alpha$, or reject it (stay at the current value) with probability $(1-\alpha)$. Repeating this thousands to tens of thousands of times makes the sample sequence converge to the posterior distribution.

Imagine it as a mechanism of randomly walking through a mountain range while staying longer in places with higher elevation (larger likelihood × prior probability).

Sharp. That's precisely the biggest bottleneck in using Bayesian calibration for CAE. The countermeasure is surrogate models, which we'll discuss in detail later, but first, let's organize the characteristics of each MCMC algorithm.

DRAM is the standard implementation in Sandia National Labs' Dakota and is the most widely used in practice. HMC/NUTS is overwhelmingly efficient if gradient information is available. For quick implementation in Python, PyMC (NUTS-based) or emcee are good choices.

This is a very important point. Using results from unconverged MCMC yields completely meaningless posterior distributions. Standard convergence diagnostics are as follows:

• $\hat{R}$ (Gelman-Rubin statistic): Variance ratio of multiple independent chains. Convergence is judged if $\hat{R} < 1.01$ (strictly $< 1.05$).
• Trace Plot: Plot the time series of the sample sequence. A "caterpillar-like" pattern with good mixing is ideal. Trends or jumps indicate non-convergence.
• Effective Sample Size (ESS): Substantial sample size considering autocorrelation. A guideline is ESS $> 400$ (for each parameter).
• Burn-in Removal: Discard initial transient samples (typically 10–50% of chain length).

The most common solution is surrogate models (proxy models). They replace the input-output relationship of FEM with a fast approximation model. Main methods are:

• Gaussian Process Regression (GP/Kriging): High affinity with the Kennedy-O'Hagan framework. Also provides prediction uncertainty.
• Polynomial Chaos Expansion (PCE): Expands the probability distribution of input parameters with orthogonal polynomials. Can be used simultaneously with uncertainty propagation.
• Neural Networks: Effective for strongly high-dimensional, nonlinear cases. However, requires large amounts of data.

For example, running FEM 100–500 times to train a GP allows each MCMC step to use only GP prediction (millisecond order). In principle, this can reduce FEM computations from $10^4 \to 10^2$.

Exactly. Therefore, surrogate validation is essential. Specifically, use Leave-one-out cross-validation to confirm prediction accuracy. $Q^2 > 0.99$ is a guideline. Another strategy is sequential surrogate updating. During MCMC, add FEM calculation points to high-probability regions of the posterior distribution and refine the surrogate. This achieves both efficiency and accuracy.