Implementing a calibration process

We provide next details on the Python implementation of a general calibration process. Later, specific examples will be provided to illustrate the concepts presented here.

As explained in this documentation, this is the list of ingredients of any calibration in ACBICI:

  1. A model to be calibrated as given by a Python function, or a text file with a rich enough sample of model evaluations.

  2. A collection of experimental points stored in a text file, with known or unknown measurement error.

  3. Prior probabilities for all the unknown parameters.

In addition, the user must decide if the covariance of the experimental error is to be estimated, if the model is to be replaced with a surrogate, and if the model discrepancy needs to be evaluated.

Defining a model

The first step in an ACBICI calibration is the definition of the model. The library defines the abstract class ACBICImodel and the user must define, for each calibration, a derived class from this one. The derived class must include two functions:

  • The initialization function (called __init__ as in any Python class) which must define the dimension of the input space and define all the unknown parameters. This value must be stored in the variable self.xdim. Then, each model parameter must be added, together with its label and a prior distribution. The syntax is, for example,

    self.addParameter(label=r'$\theta$', prior=Uniform(a=1.0, b=1.5))
    self.addParameter(label=r'$\beta$', prior=Gamma(alpha=5.0, beta=5.0))
    

    In the label field, any standard string or r-string with LaTeX symbols can be employed. Since these labels will be used in the legends of the plots, the extra burden of using LaTeX notation may pay off. Note that when defining the symbolicModel (see below), the parameters cannot be retrieved by their name. Rather, each of them is associated with the position in which it is added to the parameter list.

    For the priors, the user can employ all the priors defined in the file priors.py. In contrast with the experimental error variance (for which ACBICI provides a default prior), the user must necessarily select a prior for the parameters. The following table resumes the current priors and their required parameters

Priors

Name

Parameter

Parameter

Cauchy

mu

sigma

Gamma

alpha

beta

HalfCauchy

mu

sigma

HalfNormal

mu

sigma

Normal

mu

sigma

Uniform

a

b

Weibull

l

k

  • The function symbolicModel that, given some value of the input variables \(x\) and the parameters \(\theta\), returns the values of the model \(g(x;\theta)\).

These two functions suffice to completely define a model that can be calibrated by ACBICI.

Defining a calibration process

In a calibration process, ACBICI employs the data and functions of an ACBICImodel to optimize the posterior probabilities of the parameters. The library has been designed so that the user must create a calibrator object that takes care of all the steps of the analysis, provided the user feeds it with data and the ACBICImodel. The process for a complete calibration is always the same and follows the steps defined next:

  1. Create a calibrator object. This object, whose name is arbitrarily provided by the user, can be of one of four types:

Calibrator types

Name

Type

classicalCalibrator

A

expensiveCalibrator

B

discrepancyCalibrator

C

KOHcalibrator

D

The calibrator object is created according to its type and the user can give a name to it that will be used to write output information to disc.

  1. The experimental error in the calibrator must be defined either as fixed, or as unknown. If it is fixed, the standard deviation of its probability distribution is set with the command

    theCalibrator.setExperimentalSTDValue(0.5)
    

where here \(\sigma=0.5\). If the experimental error is unknown, a prior probability function has to be given to its standard deviation. For example, the user can issue

theCalibrator.setExperimentalSTDPrior(Gamma(alpha=2,beta=2))

and set the prior distribution to be of Gamma type. Any other prior distribution can be employed.

  1. The experimental data, which are stored in a file with the appropriate format, has to be fed into the calibrator. For example, if the data are in a file called experiments.dat, the user can load then with the commands

    experiments = np.loadtxt("data/experiments.dat")
    theCalibrator.storeExperimentalData(experiments)
    

In the case of the expensiveCalibrator and the KOHcalibrator, the analytical expression of the model is not employed (because the user has decided that it is too expensive to run). Instead, the model is sampled on pairs of input and parameter values that is deemed to cover well enough its domain. These samples are stored in a file called synthetic.dat and loaded into the calibrator with the command

sd = np.loadtxt("data/synthetic.dat")
theCalibrator.storeSyntheticData(sd)
  1. Calibrate the model. This is perform issuing the command calibrate simply as in

    theCalibrator.calibrate()
    

The calibrate command runs the algorithms corresponding to each calibration type. It has a few options that can be tuned to the particular needs:

Options for command calibrate

Name

Type

Explanation

Default value

method

string

solver method

“mcmc”

nsteps

Integer

Number of MCMC steps

10000

burn

Float

Burn-in ratio in MCMC

0.1

thin

Integer

Thinning in MCMC

1

nwalkers

Integer

Walkers in MCMC

16

At this moment, two methods can be used to solve the calibration problem, irrespective of the calibration type. The first (and default) method is “mcmc”, which stands for Markov chain Monte Carlo. This is the standard workhorse of Bayesian inference, capable to finding the posterior probability of the calibrated model by smartly sampling the parameter and hyperparameter space. In ACBICI, we use the library [emcee](https://emcee.readthedocs.io/en/stable/).

The second method available is “vbmc”, which stands for “Variational Bayesian Monte Carlo”. This is an approximate method for calculating the posterior probability distribution, picking the best one among a large, yet limited family of functions. In ACBICI, we use the library [PyVBMC](https://github.com/acerbilab/pyvbmc), which has shown to be a robust alternative to “mcmc” with remarkable time savings.

The calibration process will give some messages on the screen regarding its progress. Once finished, the user can print information on log files and obtain figures. These actions are executed with the commands

theCalibrator.printReport()
theCalibrator.plot()

The first command will write a text file in a directory with the name of the calibrator. The second command will produce some plots and has the following user options:

Options for command plot

Name

Type

Explanation

Default value

onscreen

Boolean

If true, plots are shown on the screen

True

trace

Boolean

If true, trace of MCMC is plot

True

corner

Boolean

Corner plots of posterior distributions

True

dumpfiles

Boolean

If true, write jpeg files with figures

True