Bayesian Uncertainty Quantification (UQ) on a 1D Diffusion Problem

This example demonstrates how to infer unknown model parameters given experimental observations. Uncertainties of the model parameters are inversely quantified through the Bayesian framework, which makes use of samplers including IndependentGaussianMH (Independent Gaussian Metropolis-Hastings), AffineInvariantStretchSampler or AffineInvariantDES (AffineInvariantDifferential).

Overview of Bayesian UQ

In general, assume a computational model , where are unknown model parameters. Original belief about the unknown parameters, described by the prior distribution , is usually subject to large uncertainties. Given experimental measurements of the computational model at specific configuration conditions , updated knowledge about the model parameters can be obtained by the Bayes' theorem:

where is the Likelihood function. Gaussian and TruncatedGaussian distributions are the most popular choices for the likelihood function. Variance of the Guassian or truncated Gaussian, which describes the sum of the model discrepancy and measurement error (referred to as term herein), can be either fixed a priori or inferred from the experimental measurements . The goal of Bayesian UQ is to infer the unknown parameters (and if needed) given experimental measurements , i.e., the posterior distributions .

The code architecture for performing Markov Chain Monte Carlo (MCMC) sampling in a parallelized fashion is presented in Figure 1. There are three steps for performing the sampling: proposal, model evaluation, and decision making. These three steps and the respective code objects are discussed below:

  1. Proposals: The proposals are made using the user-defined specific MCMC Samplers object which derives from the MCMC base class PMCMCBase. The specific MCMC sampler can be any variant of either serial or parallel MCMC techniques. In total, the MCMC sampler object proposes proposals which are subsequently sent to the MultiApps system.

  2. Model evaluation: The MultiApps system executes the models for each of the parallel proposals. In total, there will be model evaluations, where is the number of experimental configurations. MultiApps is responsible for distributing these jobs across a given set of processors.

  3. Decision making: The model outputs from the MultiApps system are gathered by a Reporters object to compute the likelihood function. Then, the acceptance probabilities are computed for the parallel proposals and accept/reject decisions are made. The information of accepted samples is sent to the MCMC sampler object to influence the next set of proposals.

The above three steps are repeated for a user-specified number of times.

Figure 1: Code architecture for parallelized MCMC sampling in MOOSE for performing Bayesian UQ (Dhulipala et al. (2023)).

Problem Statement

This tutorial demonstrates the application of Bayesian UQ on a 1D transient diffusion problem with Dirichlet boundary conditions on each end of the domain:

where defines the solution field. The experiments measure the average solution field across the entire domain length at the end of time . Experimental measurements are pre-generated at eight configuration points . A Gaussian noise with standard deviation of 0.05 is added to simulate the measurement error. The goal is to infer the lefthandside and righthandside Dirichlet boundary conditions, i.e., and , referred to as and for conciseness. It is worth pointing out that the term, which is the variance of the sum of the experimental measurement and model deviations, can also be inferred during the Bayesian UQ process. This tutorial will first demonstrate the inference of true boundary values by assuming a fixed variance term for the experimental error, followed by the inference of the boundary values together with variance term.

Inferring Calibration Parameters Only (Fixed Term)

Sub-Application Input

The problem defined above, with respect to the MultiApps system, is a sub-application. The complete input file for the problem is provided in Listing 1.

Listing 1: Complete input file for executing the transient diffusion problem.

left_bc = 0.13508909593042528
right_bc = -1.5530467809139854
mesh1 = 1

param1 = '${fparse left_bc}'
param2 = '${fparse right_bc}'
param3 = '${fparse mesh1}'

[Mesh<<<{"href": "../../../syntax/Mesh/index.html"}>>>]
  type = GeneratedMesh
  dim = 2
  xmax = ${param3}
  xmin = 0
  ymax = 1
  ymin = 0
  nx = 10
  ny = 10
[]

[Variables<<<{"href": "../../../syntax/Variables/index.html"}>>>]
  [u]
  []
[]

[Kernels<<<{"href": "../../../syntax/Kernels/index.html"}>>>]
  [diff]
    type = Diffusion<<<{"description": "The Laplacian operator ($-\\nabla \\cdot \\nabla u$), with the weak form of $(\\nabla \\phi_i, \\nabla u_h)$.", "href": "../../../source/kernels/Diffusion.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
  []
  [time]
    type = TimeDerivative<<<{"description": "The time derivative operator with the weak form of $(\\psi_i, \\frac{\\partial u_h}{\\partial t})$.", "href": "../../../source/kernels/TimeDerivative.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
  []
[]

[BCs<<<{"href": "../../../syntax/BCs/index.html"}>>>]
  [left]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = left
    value<<<{"description": "Value of the BC"}>>> = ${param1} # Actual = 0.15
  []
  [right]
    type = DirichletBC<<<{"description": "Imposes the essential boundary condition $u=g$, where $g$ is a constant, controllable value.", "href": "../../../source/bcs/DirichletBC.html"}>>>
    variable<<<{"description": "The name of the variable that this residual object operates on"}>>> = u
    boundary<<<{"description": "The list of boundary IDs from the mesh where this object applies"}>>> = right
    value<<<{"description": "Value of the BC"}>>> = ${param2} # Actual = -1.5
  []
[]

[Postprocessors<<<{"href": "../../../syntax/Postprocessors/index.html"}>>>]
  [average]
    type = ElementAverageValue<<<{"description": "Computes the volumetric average of a variable", "href": "../../../source/postprocessors/ElementAverageValue.html"}>>>
    variable<<<{"description": "The name of the variable that this object operates on"}>>> = u
  []
[]

[Executioner<<<{"href": "../../../syntax/Executioner/index.html"}>>>]
  type = Transient
  num_steps = 5
  dt = 0.1
  solve_type = PJFNK
  petsc_options_iname = '-pc_type -pc_hypre_type'
  petsc_options_value = 'hypre boomeramg'
[]

[Outputs<<<{"href": "../../../syntax/Outputs/index.html"}>>>]
  console<<<{"description": "Output the results using the default settings for Console output"}>>> = 'false'
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/sub.i)

Main Input

The main application, with respect to the MultiApps system, is the driver of the stochastic simulations. The main input does not perform a solve itself, couples with the sub-application for the full stochastic analysis. Details of each individual input section will is explained in details as follows.

  1. The StochasticTools block is defined at the beginning of the main application to automatically create the necessary objects for running a stochastic analysis. 

[StochasticTools<<<{"href": "../../../syntax/StochasticTools/index.html"}>>>]
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des.i)

  1. The Distributions block is used to define the prior distributions for the two stochastic boundary conditions. In this example, the two stochastic boundary conditions are assumed to follow normal distributions 

[Distributions<<<{"href": "../../../syntax/Distributions/index.html"}>>>]
  [left]
    type = Normal<<<{"description": "Normal distribution", "href": "../../../source/distributions/Normal.html"}>>>
    mean<<<{"description": "Mean (or expectation) of the distribution."}>>> = 0.0
    standard_deviation<<<{"description": "Standard deviation of the distribution "}>>> = 1.0
  []
  [right]
    type = Normal<<<{"description": "Normal distribution", "href": "../../../source/distributions/Normal.html"}>>>
    mean<<<{"description": "Mean (or expectation) of the distribution."}>>> = 0.0
    standard_deviation<<<{"description": "Standard deviation of the distribution "}>>> = 1.0
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des.i)

  1. A Likelihoods object is created to specify the type of likelihood function used in the analysis.

  • A Gaussian-type log-likelihood function is used in the current analysis, as specified by type = Gaussian and log_likelihood = true.

  • The noise term defines the experimental noise plus model deviations against experiments. Scale of the Gaussian distribution is fixed to a constant value of 0.05 through the noise_specified term, as will be explained next in the Reporters block. Note that this term can be inferred

  • The file_name argument takes in a CSV file that contains experimental values. In this example, the exp_0.05.csv file contains 8 independent experimental measurements, which are the sum of pre-simulated solution and a Gaussian noise with standard deviation of 0.05. 

[Likelihood<<<{"href": "../../../syntax/Likelihood/index.html"}>>>]
  [gaussian]
    type = Gaussian<<<{"description": "Gaussian likelihood function evaluating the model goodness against experiments.", "href": "../../../source/likelihoods/Gaussian.html"}>>>
    noise<<<{"description": "Experimental noise plus model deviations against experiments."}>>> = 'noise_specified/noise_specified'
    file_name<<<{"description": "Name of the CSV file with experimental values."}>>> = 'exp_0_05.csv'
    log_likelihood<<<{"description": "Compute log-likelihood or likelihood."}>>> = true
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des.i)

  1. In the Samplers block, AffineInvariantDES is used to sample the left and right boundary conditions according to the parallelized Markov chain Monte Carlo method at certain configuration conditions.

  • prior_distributions takes in the pre-defined prior distributions of the parameters to be calibrated.

  • num_parallel_proposals specifies the number of proposals to make and corresponding sub-applications executed in parallel.

  • initial_values defines the starting values of the parameters to be calibrated.

  • file_name takes in a CSV file that contains the configuration values. In this example, the confg.csv file contains 8 configuration values . 

[Samplers<<<{"href": "../../../syntax/Samplers/index.html"}>>>]
  [sample]
    type = AffineInvariantDES<<<{"description": "Perform Affine Invariant Ensemble MCMC with differential sampler.", "href": "../../../source/samplers/AffineInvariantDES.html"}>>>
    prior_distributions<<<{"description": "The prior distributions of the parameters to be calibrated."}>>> = 'left right'
    num_parallel_proposals<<<{"description": "Number of proposals to make and corresponding subApps executed in parallel."}>>> = 5
    file_name<<<{"description": "Name of the CSV file with configuration values."}>>> = 'confg.csv'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = PRE_MULTIAPP_SETUP
    seed<<<{"description": "Random number generator initial seed"}>>> = 2547
    initial_values<<<{"description": "The starting values of the inputs to be calibrated."}>>> = '0.1 0.1'
    previous_state<<<{"description": "Reporter value with the previous state of all the walkers."}>>> = 'mcmc_reporter/inputs'
    previous_state_var<<<{"description": "Reporter value with the previous state of all the walkers for variance."}>>> = 'mcmc_reporter/variance'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des.i)

  1. In the MultiApps block, a SamplerFullSolveMultiApp object is defined to create and run a sub-application for each sample provided by the sampler object. More specifically, a full-solve type sub-application is executed for each row of the Sampler matrix.

[MultiApps<<<{"href": "../../../syntax/MultiApps/index.html"}>>>]
  [sub]
    type = SamplerFullSolveMultiApp<<<{"description": "Creates a full-solve type sub-application for each row of each Sampler matrix.", "href": "../../../source/multiapps/SamplerFullSolveMultiApp.html"}>>>
    input_files<<<{"description": "The input file for each App.  If this parameter only contains one input file it will be used for all of the Apps.  When using 'positions_from_file' it is also admissable to provide one input_file per file."}>>> = sub.i
    sampler<<<{"description": "The Sampler object to utilize for creating MultiApps."}>>> = sample
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des.i)

  1. Under the Transfers block, SamplerReporterTransfer is utilized to transfer data from Reporters on the sub-application to a StochasticReporter on the main application. In this example, the value of a postprocessor named average, as defined in the sub-application, is transferred to a StochasticReporter named constant in the main application.

[Transfers<<<{"href": "../../../syntax/Transfers/index.html"}>>>]
  [reporter_transfer]
    type = SamplerReporterTransfer<<<{"description": "Transfers data from Reporters on the sub-application to a StochasticReporter on the main application.", "href": "../../../source/transfers/SamplerReporterTransfer.html"}>>>
    from_reporter<<<{"description": "The name(s) of the Reporter(s) on the sub-app to transfer from."}>>> = 'average/value'
    stochastic_reporter<<<{"description": "The name of the StochasticReporter object to transfer values to."}>>> = 'constant'
    from_multi_app<<<{"description": "The name of the MultiApp to receive data from"}>>> = sub
    sampler<<<{"description": "A the Sampler object that Transfer is associated.."}>>> = sample
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des.i)

  1. The Controls block makes use of MultiAppSamplerControl to pass command line arguments from the main application to the sub-applications. In this example, three arguments are passed from the main application to the sub-applications, i.e.,

  • left_bc: the Dirichlet boundary condition at the lefthandside of the domain.

  • right_bc: the Dirichlet boundary condition at the righthandside of the domain.

  • mesh1: the domain length, which is controlled by the configuration parameters in the Samplers block.

[Controls<<<{"href": "../../../syntax/Controls/index.html"}>>>]
  [cmdline]
    type = MultiAppSamplerControl<<<{"description": "Control for modifying the command line arguments of MultiApps.", "href": "../../../source/controls/MultiAppSamplerControl.html"}>>>
    multi_app<<<{"description": "The name of the MultiApp to control."}>>> = sub
    sampler<<<{"description": "The Sampler object to utilize for altering the command line options of the MultiApp."}>>> = sample
    param_names<<<{"description": "The names of the command line parameters to set via the sampled data."}>>> = 'left_bc right_bc mesh1'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des.i)

  1. Then, the Reporters section defines Reporter objects for data transfer in the MultiApps system.

  • The constant reporter of type StochasticReporter serves as a storage container for stochastic simulation results coming from the sub-application.

  • noise_specified makes use of the ConstantReporter reporter to generate a constant value, which is used for the variance term in the experimental error in this example. The value is specified to be 0.05 as defined in real_values.

  • mcmc_reporter is of type AffineInvariantDES which makes the decision whether or not to accept a proposal generated by the parallel Markov chain Monte Carlo algorithms. It takes in a likelihood function likelihoods, a sampler, and the model output from the sub-application output_value. In this example, output_value = constant/reporter_transfer:average:value means that the model output from the sub-application is obtained by a reporter named constant, which interacts with the reporter reporter_transfer to transfer the average solution field.

[Reporters<<<{"href": "../../../syntax/Reporters/index.html"}>>>]
  [constant]
    type = StochasticReporter<<<{"description": "Storage container for stochastic simulation results coming from Reporters.", "href": "../../../source/reporters/StochasticReporter.html"}>>>
  []
  [noise_specified]
    type = ConstantReporter<<<{"description": "Reporter with constant values to be accessed by other objects, can be modified using transfers.", "href": "../../../source/reporters/ConstantReporter.html"}>>>
    real_names<<<{"description": "Names for each real value."}>>> = 'noise_specified'
    real_values<<<{"description": "Values for reals."}>>> = '0.05'
  []
  [mcmc_reporter]
    type = AffineInvariantDifferentialDecision<<<{"description": "Perform decision making for Affine Invariant differential MCMC.", "href": "../../../source/reporters/AffineInvariantDifferentialDecision.html"}>>>
    output_value<<<{"description": "Value of the model output from the SubApp."}>>> = constant/reporter_transfer:average:value
    sampler<<<{"description": "The sampler object."}>>> = sample
    likelihoods<<<{"description": "Names of likelihoods."}>>> = 'gaussian'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des.i)

Available Sampling Algorithms

Currently, multiple parallel MCMC samplers have been implemented in stochastic tools module, including IndependentGaussianMH, AffineInvariantStretchSampler, and AffineInvariantDES. The above example is demonstrated on AffineInvariantDES. Switching to a different MCMC sampler only requires minor changes in the Samplers block and the Reporters block in the main input file, while the input syntax for the sub-application remains the same. Three available MCMC samplers are explained below respectively:

  • IndependentGaussianMH(IMH): Performs Metropolis-Hastings MCMC sampling with independent Gaussian proposals. Under the Samplers block, std_prop specifies the standard deviations of the independent Gaussian proposals for making the next proposal. seed_inputs takes in a reporter named mcmc_reporter which reports seed inputs values for the next proposals. The mcmc_reporter is defined under the Reporters block of type IndependentMHDecision, which performs decision making for independent Metropolis-Hastings MCMC (Calderhead, 2014).

[Samplers<<<{"href": "../../../syntax/Samplers/index.html"}>>>]
  [sample]
    type = IndependentGaussianMH<<<{"description": "Perform M-H MCMC sampling with independent Gaussian propoposals.", "href": "../../../source/samplers/IndependentGaussianMH.html"}>>>
    prior_distributions<<<{"description": "The prior distributions of the parameters to be calibrated."}>>> = 'left right'
    # previous_state = 'mcmc_reporter/inputs'
    num_parallel_proposals<<<{"description": "Number of proposals to make and corresponding subApps executed in parallel."}>>> = 5
    file_name<<<{"description": "Name of the CSV file with configuration values."}>>> = 'confg.csv'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = PRE_MULTIAPP_SETUP
    seed<<<{"description": "Random number generator initial seed"}>>> = 2547
    std_prop<<<{"description": "Standard deviations for making the next proposal."}>>> = '0.05 0.05'
    initial_values<<<{"description": "The starting values of the inputs to be calibrated."}>>> = '0.1 0.1'
    seed_inputs<<<{"description": "Reporter with seed inputs values for the next proposals."}>>> = 'mcmc_reporter/seed_input'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_imh.i)
[Reporters<<<{"href": "../../../syntax/Reporters/index.html"}>>>]
  [constant]
    type = StochasticReporter<<<{"description": "Storage container for stochastic simulation results coming from Reporters.", "href": "../../../source/reporters/StochasticReporter.html"}>>>
  []
  [noise_specified]
    type = ConstantReporter<<<{"description": "Reporter with constant values to be accessed by other objects, can be modified using transfers.", "href": "../../../source/reporters/ConstantReporter.html"}>>>
    real_names<<<{"description": "Names for each real value."}>>> = 'noise_specified'
    real_values<<<{"description": "Values for reals."}>>> = '0.05'
  []
  [mcmc_reporter]
    type = IndependentMHDecision<<<{"description": "Perform decision making for independent Metropolis-Hastings MCMC.", "href": "../../../source/reporters/IndependentMHDecision.html"}>>>
    output_value<<<{"description": "Value of the model output from the SubApp."}>>> = constant/reporter_transfer:average:value
    sampler<<<{"description": "The sampler object."}>>> = sample
    likelihoods<<<{"description": "Names of likelihoods."}>>> = 'gaussian'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_imh.i)
[Samplers<<<{"href": "../../../syntax/Samplers/index.html"}>>>]
  [sample]
    type = AffineInvariantStretchSampler<<<{"description": "Perform Affine Invariant Ensemble MCMC with stretch sampler.", "href": "../../../source/samplers/AffineInvariantStretchSampler.html"}>>>
    prior_distributions<<<{"description": "The prior distributions of the parameters to be calibrated."}>>> = 'left right'
    num_parallel_proposals<<<{"description": "Number of proposals to make and corresponding subApps executed in parallel."}>>> = 5
    file_name<<<{"description": "Name of the CSV file with configuration values."}>>> = 'confg.csv'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = PRE_MULTIAPP_SETUP
    seed<<<{"description": "Random number generator initial seed"}>>> = 2547
    initial_values<<<{"description": "The starting values of the inputs to be calibrated."}>>> = '0.1 0.1'
    previous_state<<<{"description": "Reporter value with the previous state of all the walkers."}>>> = 'mcmc_reporter/inputs'
    previous_state_var<<<{"description": "Reporter value with the previous state of all the walkers for variance."}>>> = 'mcmc_reporter/variance'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_ss.i)
[Reporters<<<{"href": "../../../syntax/Reporters/index.html"}>>>]
  [constant]
    type = StochasticReporter<<<{"description": "Storage container for stochastic simulation results coming from Reporters.", "href": "../../../source/reporters/StochasticReporter.html"}>>>
  []
  [noise_specified]
    type = ConstantReporter<<<{"description": "Reporter with constant values to be accessed by other objects, can be modified using transfers.", "href": "../../../source/reporters/ConstantReporter.html"}>>>
    real_names<<<{"description": "Names for each real value."}>>> = 'noise_specified'
    real_values<<<{"description": "Values for reals."}>>> = '0.05'
  []
  [mcmc_reporter]
    type = AffineInvariantStretchDecision<<<{"description": "Perform decision making for Affine Invariant stretch MCMC.", "href": "../../../source/reporters/AffineInvariantStretchDecision.html"}>>>
    output_value<<<{"description": "Value of the model output from the SubApp."}>>> = constant/reporter_transfer:average:value
    sampler<<<{"description": "The sampler object."}>>> = sample
    likelihoods<<<{"description": "Names of likelihoods."}>>> = 'gaussian'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_ss.i)
[Samplers<<<{"href": "../../../syntax/Samplers/index.html"}>>>]
  [sample]
    type = AffineInvariantDES<<<{"description": "Perform Affine Invariant Ensemble MCMC with differential sampler.", "href": "../../../source/samplers/AffineInvariantDES.html"}>>>
    prior_distributions<<<{"description": "The prior distributions of the parameters to be calibrated."}>>> = 'left right'
    num_parallel_proposals<<<{"description": "Number of proposals to make and corresponding subApps executed in parallel."}>>> = 5
    file_name<<<{"description": "Name of the CSV file with configuration values."}>>> = 'confg.csv'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = PRE_MULTIAPP_SETUP
    seed<<<{"description": "Random number generator initial seed"}>>> = 2547
    initial_values<<<{"description": "The starting values of the inputs to be calibrated."}>>> = '0.1 0.1'
    previous_state<<<{"description": "Reporter value with the previous state of all the walkers."}>>> = 'mcmc_reporter/inputs'
    previous_state_var<<<{"description": "Reporter value with the previous state of all the walkers for variance."}>>> = 'mcmc_reporter/variance'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des.i)
[Reporters<<<{"href": "../../../syntax/Reporters/index.html"}>>>]
  [constant]
    type = StochasticReporter<<<{"description": "Storage container for stochastic simulation results coming from Reporters.", "href": "../../../source/reporters/StochasticReporter.html"}>>>
  []
  [noise_specified]
    type = ConstantReporter<<<{"description": "Reporter with constant values to be accessed by other objects, can be modified using transfers.", "href": "../../../source/reporters/ConstantReporter.html"}>>>
    real_names<<<{"description": "Names for each real value."}>>> = 'noise_specified'
    real_values<<<{"description": "Values for reals."}>>> = '0.05'
  []
  [mcmc_reporter]
    type = AffineInvariantDifferentialDecision<<<{"description": "Perform decision making for Affine Invariant differential MCMC.", "href": "../../../source/reporters/AffineInvariantDifferentialDecision.html"}>>>
    output_value<<<{"description": "Value of the model output from the SubApp."}>>> = constant/reporter_transfer:average:value
    sampler<<<{"description": "The sampler object."}>>> = sample
    likelihoods<<<{"description": "Names of likelihoods."}>>> = 'gaussian'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des.i)

Inferring Both Calibration Parameters and the Variance Term

In addition to inferring the model parameters, the variance term can also be inferred at the same time. Both the differential evolution sampler AffineInvariantDES or stretch sampler AffineInvariantStretchSampler can be used for this purpose. Changes in the following blocks are necessary in the main input file.

  • A prior distribution need be defined under the Distributions block. In this example, a uniform distribution with wide uncertainty range is specified for the variance term as a conservative assumption.

[Distributions<<<{"href": "../../../syntax/Distributions/index.html"}>>>]
  [left]
    type = Normal<<<{"description": "Normal distribution", "href": "../../../source/distributions/Normal.html"}>>>
    mean<<<{"description": "Mean (or expectation) of the distribution."}>>> = 0.0
    standard_deviation<<<{"description": "Standard deviation of the distribution "}>>> = 1.0
  []
  [right]
    type = Normal<<<{"description": "Normal distribution", "href": "../../../source/distributions/Normal.html"}>>>
    mean<<<{"description": "Mean (or expectation) of the distribution."}>>> = 0.0
    standard_deviation<<<{"description": "Standard deviation of the distribution "}>>> = 1.0
  []
  [variance]
    type = Uniform<<<{"description": "Continuous uniform distribution.", "href": "../../../source/distributions/Uniform.html"}>>>
    lower_bound<<<{"description": "Distribution lower bound"}>>> = 0.0
    upper_bound<<<{"description": "Distribution upper bound"}>>> = 0.5
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des_var.i)
[Likelihood<<<{"href": "../../../syntax/Likelihood/index.html"}>>>]
  [gaussian]
    type = Gaussian<<<{"description": "Gaussian likelihood function evaluating the model goodness against experiments.", "href": "../../../source/likelihoods/Gaussian.html"}>>>
    noise<<<{"description": "Experimental noise plus model deviations against experiments."}>>> = 'mcmc_reporter/noise'
    file_name<<<{"description": "Name of the CSV file with experimental values."}>>> = 'exp_0_05.csv'
    log_likelihood<<<{"description": "Compute log-likelihood or likelihood."}>>> = true
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des_var.i)
[Reporters<<<{"href": "../../../syntax/Reporters/index.html"}>>>]
  [constant]
    type = StochasticReporter<<<{"description": "Storage container for stochastic simulation results coming from Reporters.", "href": "../../../source/reporters/StochasticReporter.html"}>>>
  []
  [mcmc_reporter]
    type = AffineInvariantDifferentialDecision<<<{"description": "Perform decision making for Affine Invariant differential MCMC.", "href": "../../../source/reporters/AffineInvariantDifferentialDecision.html"}>>>
    output_value<<<{"description": "Value of the model output from the SubApp."}>>> = constant/reporter_transfer:average:value
    sampler<<<{"description": "The sampler object."}>>> = sample
    likelihoods<<<{"description": "Names of likelihoods."}>>> = 'gaussian'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des_var.i)

  • Finally, in the Samplers block, the AffineInvariantDES sampler takes in the prior distribution of the variance term via the prior_variance argument. In this example, the prior distribution of the variance term has been defined in variance under the Distributions block.

[Samplers<<<{"href": "../../../syntax/Samplers/index.html"}>>>]
  [sample]
    type = AffineInvariantDES<<<{"description": "Perform Affine Invariant Ensemble MCMC with differential sampler.", "href": "../../../source/samplers/AffineInvariantDES.html"}>>>
    prior_distributions<<<{"description": "The prior distributions of the parameters to be calibrated."}>>> = 'left right'
    num_parallel_proposals<<<{"description": "Number of proposals to make and corresponding subApps executed in parallel."}>>> = 5
    file_name<<<{"description": "Name of the CSV file with configuration values."}>>> = 'confg.csv'
    execute_on<<<{"description": "The list of flag(s) indicating when this object should be executed. For a description of each flag, see https://mooseframework.inl.gov/source/interfaces/SetupInterface.html."}>>> = PRE_MULTIAPP_SETUP
    seed<<<{"description": "Random number generator initial seed"}>>> = 2547
    initial_values<<<{"description": "The starting values of the inputs to be calibrated."}>>> = '0.1 0.1'
    previous_state<<<{"description": "Reporter value with the previous state of all the walkers."}>>> = 'mcmc_reporter/inputs'
    previous_state_var<<<{"description": "Reporter value with the previous state of all the walkers for variance."}>>> = 'mcmc_reporter/variance'
    prior_variance<<<{"description": "The prior distribution of the variance parameter to be calibrated."}>>> = 'variance'
  []
[]
(moose/modules/stochastic_tools/test/tests/samplers/mcmc/main_des_var.i)

Result and Analysis

The output JSON file contains samples across different processors at different timesteps. Practically, all MCMC samples are generated with serial auto-correlations. Therefore, it is essential to diagnose the sample auto-correlation for predictable UQ quality. Two primary diagnostics are presented here for analyzing the sample quality: the integrated auto-correlation time (), and the effective sample size (ESS) Dhulipala et al. (2023). describes the interval after which a sample can be considered independent given a current sample index, and ESS describes the number of "effective" samples after taking into consideration the sample auto-correlation. To obtain "independent" samples from the Markov chains, the first 100 samples are discarded for "burn-in", with the remaining samples thinned by . The ESS of the generated samples are calculated based on the thinned samples. As presented in Table 1, for IMH, SS and DES decreases, and ESS for the three algorithms increases correspondingly, indicating the best sampling qualify by DES, followed by SS, then IMH.

Table 1: Integrated autocorrelation time () and effective sample size (ESS) when inferring model parameters

AlgorithmESS
IMH123.9089.528
SS41.77946.719
DES11.217676.547

When inferring only the boundary conditions, all three samplers are applied, with 1,000 serial steps simulated with 5 parallel proposals per step. Figure 2 shows the posterior distributions of the calibration parameters and . A strong linearity is observed for the two parameters as would be expected from the underlying physics. The IMH samples appear to be extremely localized compared to DES. The SS samples are almost as good as samples generated by DES in terms of exploring the parameter space. The DES samples are the best from this aspect. More detailed posterior distributions are presented in Figure 3 for all three algorithms IMH, SS and DES. It is generally observed that the posteriors using IMH are restricted in the parameter space. Whereas, the posteriors for both SS and DES are spread across the parameter space and look comparable to each other. This can be observed by comparing Figure 3b and Figure 3c for the posterior of , and Figure 3e and Figure 3f for the posterior of .

Figure 2: Scatter plot of the samples from IMH, SS, and DES samplers when only inferring the model parameters.

Figure 3: Posterior distributions for the diffusion time derivative problem when inferring only the model parameters. (a), (b) and (c) are the posteriors using IMH, SS, and DES samplers, respectively. (d), (e) and (f) are the posteriors using IMH, SS, and DES samplers, respectively.

When inferring both boundary conditions and the term, only the SS and DES samplers are considered given the poor performance of the IMH sampler evidenced from the above case. Again, 1,000 serial steps are simulated with 5 parallel proposals per step. Sampling diagnostics are presented for the two samplers in Table 2. The for both samplers are comparable, and the for DES increased compared to the previous case of inferring only model parameters, possible caused by the increased problem complexity due to the inclusion of the term. For the ESS, the DES shows a three time enhancement compared to SS, due to the higher residual sample autocorrelation upon chain thinning.

Table 2: Integrated autocorrelation time () and effective sample size (ESS) when inferring both model parameters and the term.

AlgorithmESS
SS36.22851.991
DES45.542155.0

Scatter plot of the inferred model parameters and is presented in Figure 4. Again, the two model parameters exhibit almost a perfect correlation which is consistent with the observation in the previous case, due to the ill-posedness of the problem.

Figure 4: Scatter plot of the samples from IMH, SS, and DES samplers when inferring both model parameters and the term.

Detailed posterior distributions of , and the term are presented in Figure 5 for the SS and DES samplers. Both samplers are producing quite consistent posterior distributions. Specifically, Figure 5a and Figure 5d need to be compared for , Figure 5b and Figure 5e for , and Figure 5c and Figure 5f for the term. Both samplers produce consistent posterior distributions, and the posterior distributions of the term are approximately centered around the true value of the Gaussian noise (denoted by the solid vertical line) that was added to replicate the “experimental” data.

Figure 5: Posterior distributions for the diffusion time derivative problem when inferring the model parameters and the term. (a), (b) and (c) show the posteriors using SS. (d), (e) and (f) show the posteriors using DES.

As shown in Figure 3 and Figure 5, posteriors of the model parameters and potentially the term are much more constrained compared to the priors, reflecting the integration of experimental information into the models through the Bayesian UQ process. Followup forward UQ based on the updated posteriors will provide more accurate information about the output uncertainty.

References

  1. C. J. T. Braak. A markov chain monte carlo version of the genetic algorithm differential evolution: easy bayesian computing for real parameter spaces. Statistics and Computing, 16:239–249, 2006. doi:10.1007/s11222-006-8769-1.[BibTeX]
  2. B. Calderhead. A general construction for parallelizing metropolis-hastings algorithms. Proceedings of the National Academy of Sciences, 111(49):17408–17413, 2014. doi:10.1073/pnas.1408184111.[BibTeX]
  3. Som LakshmiNarasimha Dhulipala, Daniel Schwen, Yifeng Che, Ryan Terrence Sweet, Aysenur Toptan, Zachary M Prince, Peter German, and Stephen R Novascone. Massively parallel bayesian model calibration and uncertainty quantification with applications to nuclear fuels and materials. Technical Report, Idaho National Laboratory (INL), Idaho Falls, ID (United States), 2023.[BibTeX]
  4. J. Goodman and J. Weare. Ensemble samplers with affine invariance. Communications in applied mathematics and computational science, 5(1):65–80, 2010. doi:10.2140/camcos.2010.5.65.[BibTeX]
  5. B. Nelson, E. B. Ford, and M. J. Payne. Run dmc: an efficient, parallel code for analyzing radial velocity observations using n-body integrations and differential evolution markov chain monte carlo. The Astrophysical Journal Supplement Series, 11:11–25, 2013. doi:10.1088/0067-0049/210/1/11.[BibTeX]
  6. M. D. Shields, D. G. Giovanis, and V. S. Sundar. Subset simulation for problems with strongly non-gaussian, highly anisotropic, and degenerate distributions. Computers & Structures, 245:106431, 2021. doi:10.1016/j.compstruc.2020.106431.[BibTeX]