Skip to content

Bayesian perspective on model calibration

Note

This section focuses on stochastic models, but the principles also apply to deterministic models with simplifications.

From a statistical perspective, we are interested in the joint distribution of model parameters and latent states after seeing the observed data. In other words, the distribution we are interested in is:

\[ p(\theta, X \mid y^{\text{obs}}). \]

We have dropped the \(0:T\) subscript for notational brevity. As typical in Bayesian analyses, we cannot sample from this distribution directly. Instead, we use Bayes' rule to write:

\[ p(\theta, X \mid y^{\text{obs}}) = \frac{ p(y^{\text{obs}} \mid X, \theta) p(X \mid \theta) p(\theta)}{ p(y^{\text{obs}}) }. \]

In the above expansion, \(p(y^{\text{obs}} \mid X, \theta)\) is the likelihood function, \(p(X \vert \theta)\) contains the model dynamics, and \(p(\theta)\) is the prior distribution over parameters. The term \(p(y^{\text{obs}})\) can be viewed as a normalizing constant that is independent from \(X\) and \(\theta\). The likelihood is a conditional probability density, but we view it as a function of \(X\) and \(\theta\) considering that \(y^{\text{obs}}\) is known,

\[ \mathcal L(X, \theta) \triangleq p(y^{\text{obs}} \mid X, \theta). \]

The likelihood quantifies how well the model output matches the observed data for a given set of parameters and latent state trajectory.

For deterministic models, the main difference is that the latent state trajectory \(X\) is fully determined by the input parameters \(\theta\). Therefore, we can write \(x = X(\theta)\) and simplify the likelihood function to depend only on \(\theta\). A variety of methods exist to optimize or sample from the posterior distribution of \(\theta\) given the data in this case.

Obtaining posterior samples for stochastic models

For agent-based models, we do not have a closed form equation describing the transition density \(p(X \mid \theta)\). However, what we can do is simulate trajectories \(\tilde x\) and outputs \(\tilde y\) given parameters \(\theta\) by running the computer model.

Sampling algorithms are often differentiated by their loss function and approach stages.

Loss: In computational epidemiology, we rarely know the true likelihood function, \(p(y \mid X, \theta)\) exactly. Therefore, we either compose a pseudo-likelihood function based on the nature of the data (e.g., normal, binomial, Poisson) or take a likelihood-free approach that avoids explicit evaluation of the likelihood function. We will discuss each of these approaches in more detail in the next section.

Stages: The goal is to obtain samples from the posterior of the joint distribution \(p(X, \theta \mid y^{\text{obs}})\). Some methods use a single stage to jointly sample pairs \((X, \theta)\) directly from the posterior. Other methods approach the problem in a two-stage conditional manner by first sampling parameters from the marginal posterior \(p(\theta \mid y^{\text{obs}})\) and then sampling latent state trajectories from the "smoothing" distribution \(p(X \mid y^{\text{obs}}, \theta)\).

Brief introduction to algorithms used in calibration

Many algorithms have been developed to obtain posterior samples. Here, we briefly introduce a few algorithms that are used frequently in epidemiological modeling applications.

  • Markov Chain Monte Carlo (MCMC): Sequentially propose parameters and accept or reject based on loss. Requires many sequential simulations, including a "burn-in" period, to ensure the accepted samples are distributed according to the target distribution. A possible speed-up, called delayed acceptance, is to train an emulator to predict the likelihood given input parameters and only run the real simulation if the acceptance probability is sufficiently high.

  • Approximate Bayesian Computation (ABC): Instead of composing a pseudo-likelihood function, use a kernel (often the \(L_2\) norm) to evaluate the distance between observed and simulated data. Keep samples that are within distance \(\varepsilon\).

  • Sampling with Importance Resampling: To sample from distribution \(p\), draw samples from an "easy to sample" distribution \(q\), weight each likelihood by the importance weight \(p/q\), and keep samples with probability proportional to the weighted likelihood.

  • Sequential Monte Carlo (SMC): Also called particle filtering. Similar to sampling with importance resampling, but offers a lower variance estimate of the likelihood in exchange for complexity. SMC^2 runs SMC in parameter space using random walk exploration while also running SMC to evaluate likelihood at each point.

  • History Matching: Iteratively remove implausible regions of parameter space using model output emulation. Can be used as a pre-processing step to other sampling methods.

  • Simulation-Based Inference: Train a neural network to learn a mapping from observed data \(y\) to parameters of a distribution describing the model input parameters. For example, the network might learn the mean and covariance parameters of a multi-variate normal distribution. Once the model has been trained, the real observation data can be fed in to reveal the parametric posterior distribution over model inputs.

  • Bayesian Synthetic Likelihood: Assumes that the likelihood can be well approximated by a parametric distribution, often a multi-variate normal. Run many simulations for each input and use results to estimate the likelihood distribution parameters.