Continuous Calibration

assess_convergence(stepsize::Array{Float64},acceptance::Array{Float64},nx::Int,ntheta::Int)

Function to assess the convergence of the M-H stepsize optimization results.

Positional arguments

• stepsize::Array{Float64} A n by nx+ntheta Array containing the M-H step sizes used in the stepsize optimization algorithm.
• acceptance::Array{Float64} A n by nx+ntheta Array containing the calculated M-H acceptance rates corresponding to stepsize.
• nx::Int The number of x dimesions.
• ntheta::Int The number of θ dimensions.

Returns

• secants::Tuple{Array{Float64}} A Tuple of length 2, with each element being a (n-1) by nx+ntheta Array.
  • The Arrays contain the calculated secants of the acceptance rates and stepsizes, respectively, calculated as a function of the epochs, relative to the final epoch.

Details This function calculates the slope of the secant line between each element of the acceptance and stepsize Arrays relative to the last element.

auto_stepsize(model,data::DataStr,nbatch::Int,batchsize::Int,prior_data::PriorData,
nx::Int,ntheta::Int,nobs::Int,theta_init::Vector{Float64},
init::Float64,target::Tuple{Float64},eta::Float64,factor::Float64,offset::Float64)

Function to calculate the appropriate step size for the Metropolis-Hastings algorithm, for a target acceptance ratio.

Positional arguments

• model Surrogate model.
• data::DataStr Struct containing the computer simulator and experimental data.
• nbatch::Int The number of batches of MCMC simulations to run.
• batchsize::Int The number of MCMC iterations ro run per batch.
• prior_data::PriorData Struct containing information on the variables' prior distributions.
• nx::Int The number of x dimensions.
• ntheta::Int The number of θ dimensions.
• nloc::Int The number of unique settings of x in the experimental data.
• theta_init::Vector{Float64} The settings of θ at which to initialize the MCMC.
• init::Float64 The inital step size for the start of this algorithm.
• target::Tuple{Float64} A length 2 Tuple containing the target acceptance rates for θ and ρ, respectively.
• scale::Float64 Scaling parameter to pass into the stepsize_adjust function.
• shape::Float64 Shape parameter to pass into the stepsize_adjust function.
• offset::Float64 Offset parameter to pass into the stepsize_adjust function.

Returns

• stepsize::StepSize Struct containing the calculated stepsizes that will result in the target acceptance rate.
• stepsize_hist::Array{Float64} An Array containing the historic data of the stepsizes used in this algorithm.
• acceptance_hist::Array{Float64} An Array containing the historic data of the acceptance rates used in this algorithm.

Details This algorithm runs nbatch batches of MCMC with M-H updates, each batch having a length of batchsize. The algorithm starts at a step size specified by init for the first batch. After each batch, the algorithm calculates the average acceptance rate of all M-H samples from all preceding batches. This value is passed to update_stepsize to calculate the proposed adjustment to the step size for the next batch. The step size for the next batch is then calculated as the average of all previous step sizes and the value calculated from update_stepsize.

find_stepsize(model,data::DataStr,nbatch::Int,batchsize::Int,prior_data::PriorData,
nx::Int,ntheta::Int,nloc::Int)
find_stepsize(model,data::DataStr,nbatch::Int,batchsize::Int,prior_data::PriorData,
nx::Int,ntheta::Int,nloc::Int;theta_init::Vector{Float64},method::Int,make_plots::Bool,
show_plots::Bool,save_plots::Bool,init::Float64,target::Tuple{Float64},scale::Float64,
shape::Float64,offset::Float64,mdl_apnd::String)

Function to solve for the appropriate step size for θ and ρ in the Metropolis-Hastings algorithm for Bayesian calibration.

Keyword arguments

• model Surrogate model
• data::DataStr Struct containing the computer simulator and experimental data.
• nbatch::Int The number of batches of MCMC simulations to run.
• batchsize::Int The number of MCMC iterations ro run per batch.
• prior_data Struct containing parameter prior distribution informaiton.
• nx::Int The number of x dimensions.
• ntheta::Int The number of θ dimensions.
• nloc::Int The number of unique settings of x in the experimental data.

Keyword arguments

• theta_init::Union{Vector{Float64},Float64} The settings of θ at which to initialize the MCMC.
  • default value of 0.5 is used for each dimension of theta if not specified.
• make_plots::Bool Indicator of whether to generate plots from the results of the algorithm.
  • default value of true
• show_plots::Bool Indicator of whether to display the plots generated.
  • default value of true
• save_plots::Bool Indicator of whether to save the plots generated.
  • default value of true
• init::Float64 The inital step size for the start of this algorithm.
  • By default, 1E-3 is used.
• target::Vector{Float64} A length 2 Tuple containing the target acceptance rates for θ and ρ, respectively.
  • By defualt, 0.3 is used for both.
• scale::Float64 Scaling parameter to pass into the stepsize_adjust function.
  • By default, 2.0 is used.
• shape::Float64 Shape parameter to pass into the stepsize_adjust function.
  • By defualt, 10.0 is used.
• offset::Float64 Offset parameter to pass into the stepsize_adjust function.
  • By default, 1.5 is used.
• mdl_apnd::String String to append to the beginning of generated plot file names.

Returns

• stepsize::StepSize Struct containing the calculated stepsizes that will result in the target acceptance rate.

plot_stepsize_opt(stepsize::Array{Float64},acceptance::Array{Float64},nx::Int,
ntheta::Int,show_plots::Bool,save_plots::Bool,mdl_apnd::String)

Function to plot the results of the M-H stepsize optimization algorithm.

Positional arguments

• stepsize::Array{Float64} A n by nx+ntheta Array containing the M-H step sizes used in the stepsize optimization algorithm.
• acceptance::Array{Float64} A n by nx+ntheta Array containing the calculated M-H acceptance rates corresponding to stepsize.
• nx::Int The number of x dimesions.
• ntheta::Int The number of θ dimensions.
• show_plots::Bool An indication of whether the plots should be saved.
• save_plots::Bool An indication of whether the plots should be displayed.
• mdl_apnd::String String to append to the front of the plots' file names.

stepsize_adjust(diff::Float64,scale::Float64,shape::Float64,offset::Float64)

Function to adjust a given stepsize based on the difference between its calculated acceptance rate and the target.

Positional arguments

• diff::Float64 The difference between the calculated acceptance ratio and the target.
• scale::Float64 Scaling factor for the equation.
• shape::Float64 Shape parameter for the equation.
• offset::Float64 Offset parameter for skewing the adjustment to correct greater in the positive or negative direction

Returns

• factor::Float64 The adjustment by which to multiply the current step size.

Details The function uses an arc-tangent-based form for calculating the adjustment. The arc-tangent is calculated for the product of the difference and the shape. The arc-tangent is them multiplied by scale This value is then shifted based on the offset

• positive values for offset correspond to greater corrections of negative differences (the acceptance ratio is too small)
• negative values for offset correspond to greater corrections of positive differences (the acceptance ratio is too large)

The exponentiation of this value is returned

Deprecated stepsizegd(data::DataStr,nbatch::Int,nsize::Int,priordata::PriorData,nx::Int,ntheta::Int,nobs::Int,thetainit::Vector{Float64}) stepsizegd(data::DataStr,nbatch::Int,nsize::Int,priordata::PriorData,nx::Int,ntheta::Int,nobs::Int,thetainit::Vector{Float64},init::Float64,target::Tuple{Float64},eta::Float64) Function to optimize the M-H step size via gradient descent.

Positional arguments

• data::DataStr Struct containing the computer simulator and experimental data.
• nbatch::Int The number of batches of MCMC simulations to run.
• batchsize::Int The number of MCMC iterations ro run per batch.
• nx::Int The number of x dimensions.
• ntheta::Int The number of θ dimensions.
• nobs::Int The number of unique settings of x in the experimental data.
• theta_init::Vector{Float64} The settings of θ at which to initialize the MCMC.

Optional arguments

• init::Float64 The inital step size for the start of this algorithm.
  • By default, 1E-3 is used.
• target::Tuple{Float64} A length 2 Tuple containing the target acceptance rates for θ and ρ, respectively.
  • By defualt, 0.3 is used for both.
• eta::Float64 Learning rate for the gradient descent.
  • By default, 0.1 is used.

Returns

• stepsize::StepSize Struct containing the calculated stepsizes that will result in the target acceptance rate.
• stepsize_hist::Array{Float64} An Array containing the historic data of the stepsizes used in this algorithm.
• acceptance_hist::Array{Float64} An Array containing the historic data of the acceptance rates used in this algorithm.

Details This algorithm runs nbatch batches of MCMC with M-H updates, each batch having a length of batchsize. The algorithm starts at a step size specified by init for the first batch. After each batch, the algorithm calculates the acceptance rate of that batch. The gradient of a SSE objective function (as a function of the batch acceptance rate compared to the target) is calculated as a function of the log step size of θ and ρ. The the log step size for each θ and ρ is updated by subtracting the gradient multiplied by eta.

update_stepsize!(acceptance::Vector{Float64},stepsize::StepSize,history::Array{Float64},
nx::Int,ntheta::Int,target::Vector{Float64},scale::Float64,shape::Float64,offset::Float64,
weight::Float64)

Function to update the Metropolis-Hastings algorithm step sizes based on the calculated acceptance ratios compared to their target values.

Positional arguments

• acceptance::Vector{Float64} Array containing the acceptance value for each M-H update.
• stepsize::StepSize Struct containing the step sizes used during the M-H updates for θ and ρ.
• history::Array{Float64} History of previous step sizes.
• ntheta::Int The number of θ variables using M-H updates.
• nx::Int The number of ρ variables using M-H updates.
• target::Vector{Float64} A length 2 Tuple containing the target acceptance rates for θ and ρ, respectively.
• scale::Float64 Scaling parameter to pass into the stepsize_adjust function.
• shape::Float64 Shape parameter to pass into the stepsize_adjust function.
• offset::Float64 Offset parameter to pass into the stepsize_adjust function.
• weight::Float64 weight to use for the weighted average.

weighted_avg(Values::Vector{Bool},new::Vector{Bool},weight::Float64)

Function to calculate the weighted average of a set of data. Implementation for weighted average of acceptance values. This function assumes a weight of unity for the old values in values and uses weight for new

Positional arguments

• values::Vector{Bool} Old acceptance values.
• new::Vector{Bool} New acceptance values.
• weight::Float64 weight for the new values.

weighted_avg(values::Vector{Float64},new::Float64,weight::Float64)

Function to calculate the weighted average of a set of data. Implementation for weighted average of step size values. This function assumes a weight of unity for the old values in values and uses weight for new

Positional arguments

• values::Vector{Float64} Old values for the average.
• new::Float64 New value added to the average.
• weight::Float64 weight for the new value.

expit(c1::Float64,c2::Float64,gamma::Float64)

Function to transform the variable of γ ∈ in [-∞,∞] to [c2,c2+c1].

Positional arguments

• c1::Float64 Scale parameter of original variable space. If the variable is defined in [A,B], this parameter is B-A.
• c2::Float64 Minimum parameter of the original variable space. If the variable is defined in [A,B], this parameter is A.
• gamma::Float64 Value to be converted to the [A,B] interval.

Returns

• x A scalar Float transformed into [A,B] using c1*exp(γ)/(1+exp(γ)-c2)

gibbs_delta(data::DataStr,sig2::Float64,tau2::Float64,eta::Vector{Float64},
corr::Array{Float64,2},nloc::Int,nrep::Int)

Function to draw a sample from the posterior distribution of the discrepancy term, δ.

Positional arguments

• data::DataStr Data structure containing the computer simulator and experimental data.
• sig2::Float64 Sample of the data model error variance.
• tau2::Float64 Sample of the discrepancy covariance matrix variance parameter.
• corr::Array{Float64,2} Correlation matrix of the discrepancy covariance.
• nloc::Int The number of locations at which the data is observed.
• nrep::Int The number of independent observations of the data at the locations from nloc.

Returns

• sample A Vector of length n containing the posterior draw of δ.

Details The full conditional distribution of δ is proportional to the multivariate normal likelihood function multiplied by the multivarite normal prior of δ. This can be simplified to a well known multivariate normal posterior distribution. Given p(δ) = MVN(0,Σ) and L(.|x,y) = ∏MVN(η+δ,σ^2I), the posterior distribution for δ can be solved as the following: p(δ|.) ∼ MVN(bnAn,An) where An = (Σ^-1 + m^2(σ^2I)^-1)^-1 bn = m^2((σ^2I)^-1)*(ȳ-η) m is the number of independent observations of the multivariate normal data.

gibbs_sig2(prior_data::PriorData,data::DataStr,eta::Float64,nrep::Int)

Function to draw a sample from the posterior distribution of the data model variance term (σ^2). Implementation for a univariate normal distribution.

Positional arguments

• prior_data::PriorData Data structure containing the prior distribution parameters.
• data::DataStr Data structure containing the computer simulator and experimental data.
• eta::Float64, Surrogate model prediction.
• nrep::Int The number of independent observations of the multivariate normal distribution data model.

Returns

• sample A scalar float containing the new sample of σ^2 from its posterior distribution.

Details The full conditional distribution of σ^2 is proportional to a normal likelihood function multiplied by an inverse gamma distribution (the σ^2 prior distribution). This can be simplified to a well known inverse gamma posterior distribution form for Gibbs updates. Given prior distribution parameters for σ^2 of IG(α,β) and likelihood function of ∏^m(N(η,σ^2)), the posterior distribution is solved as the following: p(τ^2|.) ∼ IG(α+m/2,β+0.5*SSE) Where m is the number of independent observations from this data model, and SSE is the sum of squared errors between the data and the surrogate model.

gibbs_sig2(prior_data::PriorData,data::DataStr,eta::Vector{Float64},delta::Vector{Float64},
lik_power::Int,nrep::Int)

Function to draw a sample from the posterior distribution of the data model variance term (σ^2). Implementation for a multivariate normal distribution.

Positional arguments

• prior_data::PriorData Data structure containing the prior distribution parameters.
• data::DataStr Data structure containing the computer simulator and experimental data.
• eta::Vector{Float64}, Vector containing surrogate model prediction.
• delta::Vector{Float64} Vector containing discrepancy function sample.
• lik_power::Int Total number of data observations (number of measurement locations times the number of independent samples of those locations).
• nrep::Int The number of independent observations of the multivariate normal distribution data model.

Returns

• sample A scalar float containing the new sample of σ^2 from its posterior distribution.

Details The full conditional distribution of σ^2 is proportional to a multivariate normal likelihood function multiplied by an inverse gamma distribution (the σ^2 prior distribution). This can be simplified to a well known inverse gamma posterior distribution form for Gibbs updates. Given prior distribution parameters for σ^2 of IG(α,β) and likelihood function of ∏MVN(η+δ,σ^2I), the posterior distribution is solved as the following: p(τ^2|.) ∼ IG(α+nm/2,β+0.5λ'λ) Where n is the length of the multivariate normal distribution in the data model, m is the number of independent observations from this data model, and λ is the vector y_i - η - δ.

gibbs_tau2(prior_data::PriorData,delta::Vector{Float64},corr::Array{Float64,2},nloc::Int64)

Function to draw a sample from the posterior distribution of the discrepancy variance term (τ^2).

Positional arguments

• prior_data::PriorData Data structure containing the prior distribution parameters.
• delta::Vector{Float64} Vector containing a sample of the discrepancy function.
• corr::Array{Float64,2} Correlation matrix of the discrepancy term.
• nobs::Int64 Integer describing the number of data points in a single independent observation of the data model's multivariate normal distribution.

Returns

• sample A scalar Float containing the new sample of τ^2 from its posterior distribution.

Details The full conditional distribution of τ^2 is proportional to a multivariate normal distribution (the discrepancy term's prior distribution) multiplied by an inverse gamma distribution (the τ^2 prior distribution). This can be simplified to a well known inverse gamma posterior distribution form for Gibbs updates. Given prior distribution parameters for τ^2 of IG(α,β) and prior distribution for δ of MVN(0,τ^2C), the posterior distribution is solved as the following: p(τ^2|.) ∼ IG(α+n/2,β+0.5δ'C^(-1)δ) Where n is the length of the discrepancy term (the same length as the data model's multivariate normal distribution).

logit(c1::Float64,c2::Float64,var::Float64)

Function to transform a variable defined in [c2,c1+c2] to γ space defined in [-∞,∞].

Positional arguments

• c1::Float64 Scale parameter of original variable space. If the variable is defined in [A,B], this parameter is B-A.
• c2::Float64 Minimum parameter of the original variable space. If the variable is defined in [A,B], this parameter is A.
• var::Float64 Variable to be converted to the [-∞,∞] interval.

Returns

• gamma A scalar Float transformed into [-∞,∞] by scaling x to [0,1] with x'=(x-c2)/c1 and then using the logit function, ln(x'/(1-x'))

metropolis_rho(prior_data::PriorData,data::DataStr,rho::Vector{Float64},delta::Vector{Float64},
tau2::Float64,k::Int64,stepsize::Float64,nx::Int64,nloc::Int64)

Function to draw a sample of ρ's posterior distribution using the Metropolis-Hastings algorithm.

Positional arguments

• prior_data::PriorData Data structure containing the prior distribution parameters.
• data::DataStr Data structure containing the computer simulator and experimental data.
• rho::Vector{Float64} Vector containit ρ, whose value at index k will be updated.
• delta::Vector{Float64} Vector of a sample of the discrepancy function.
• tau2::Float64 Discrepancy covariance matrix variance parameter.
• k::Int64 Indexing integer prescribing which ρ dimension is being sampled.
• stepsize::Float64 Prescribed stepsize to use in the proposal distribution.
• nx::Int64 Integer describing the number of independent control variables (x).
• nloc::Int64 Integer describing the number of data points in a single independent observation of the data model's multivariate normal distribution.

Returns

• output::MetropolisInfo A data struct containing the value sample from the Metropolis-Hastings algorithm, the acceptance probability for the proposed value, and the Boolean value indicating if the proposed value was accepted or rejected.
• corr::Array{Float64,2} Calculated correlation matrix for the accepted value of ρ.

Details The current value of ρ[k] is transformed to the real number line using the expit function. A new value is proposed on the real number line from a normal distribution using this transformed value as the mean and a standard deviation equal to stepsize. The proposed value is transformed back to the ρ space, ρ. The proportional posterior values are calculated for the current value of ρ and the proposed value, ρ. These are denoted as π(ρ) and π(ρ). The proposal probabilities are calculated using ρ and ρ, denoted as J(ρ) and J(ρ). The acceptance probability of ρ is calculated as π(ρ)/π(ρ)J(ρ*)/J(ρ). A random sample form U(0,1) determines if the proposed value is accepted or rejected.

metropolis_theta(model,prior_data::PriorData,data::DataStr,theta::Vector{Float64},
sig2::Float64,k::Int64,stepsize::Float64)

Function to draw a sample of θ's posterior distribution using the Metropolis-Hastings algorithm. This implementation is for a univariate normal distribution data model.

Positional arguments

• model Surrogate model for predicting response.
• prior_data::PriorData Data structure containing the prior distribution parameters.
• data::DataStr Data structure containing the computer simulator and experimental data.
• theta::Vector{Float64} Vector of θ values, whose value at index k will be updated.
• sig2::Float64 Data model variance parameter, used for calculating the likelihood.
• k::Int64 Indexing integer prescribing which θ dimension is being sampled.
• stepsize::Float64 Prescribed stepsize to use in the proposal distribution.

Returns

• output::MetropolisInfo A data struct containing the value sample from the Metropolis-Hastings algorithm, the acceptance probability for the proposed value, and the Boolean value indicating if the proposed value was accepted or rejected.
• eta::Vector{Float64} Vector of surrogate model prediction for the accepted θ value

Details The current value of θ[k] is transformed to the real number line using the expit function. A new value is proposed on the real number line from a normal distribution using this transformed value as the mean and a standard deviation equal to stepsize. The proposed value is transformed back to the θ space, θ. The proportional posterior probabilities are calculated for the current value of θ and the proposed value, θ. These are denoted as π(θ) and π(θ). The proposal probabilities are calculated using θ and θ, denoted as J(θ) and J(θ). The acceptance probability of θ is calculated as π(θ)/π(θ)J(θ*)/J(θ). A random sample form U(0,1) determines if the proposed value is accepted or rejected.

metropolis_theta(model,prior_data::PriorData,data::DataStr,theta::Vector{Float64},
delta::Vector{Float64},sig2::Float64,k::Int64,stepsize::Float64)

Function to draw a sample of θ's posterior distribution using the Metropolis-Hastings algorithm. This implementation is for a multivariate normal distribution data model.

Positional arguments

• model Surrogate model for predicting response.
• prior_data::PriorData Data structure containing the prior distribution parameters.
• data::DataStr Data structure containing the computer simulator and experimental data.
• theta::Vector{Float64} Vector of θ values, whose value at index k will be updated.
• delta::Vector{Float64} Vector of discrepancy values, used for calculating the likelihood.
• sig2::Float64 Data model variance parameter, used for calculating the likelihood.
• k::Int64 Indexing integer prescribing which θ dimension is being sampled.
• stepsize::Float64 Prescribed stepsize to use in the proposal distribution.

Returns

• output::MetropolisInfo A data struct containing the value sample from the Metropolis-Hastings algorithm, the acceptance probability for the proposed value, and the Boolean value indicating if the proposed value was accepted or rejected.
• eta::Vector{Float64} Vector of surrogate model prediction for the accepted θ value

Details The current value of θ[k] is transformed to the real number line using the expit function. A new value is proposed on the real number line from a normal distribution using this transformed value as the mean and a standard deviation equal to stepsize. The proposed value is transformed back to the θ space, θ. The proportional posterior probabilities are calculated for the current value of θ and the proposed value, θ. These are denoted as π(θ) and π(θ). The proposal probabilities are calculated using θ and θ, denoted as J(θ) and J(θ). The acceptance probability of θ is calculated as π(θ)/π(θ)J(θ*)/J(θ). A random sample form U(0,1) determines if the proposed value is accepted or rejected.

mcmc!(model,data::DataStr,prior_data::PriorData,sample_vals::BulkVarsStruct,start::Int,
stop::Int,stepsize::StepSize,nx::Int,ntheta::Int,nloc::Int)

Function to perform the Markov chain Monte Carlo simulation to draw from the posterior distributions of the parameters.

Positional arguments

• model Surrogate model of computer simulator.
• data::DataStr Struct containing the computer simulator and experimental data.
• prior_data::PriorData Struct containing the prior distribution data for the model parameters.
• sample_vals::BulkVarsStruct The struct to store the sampling results.
• start::Int The index in sample_vals at which to start the sampling.
• stop::Int The index in sample_vals at which to stop the sampling.
• stepsize::StepSize Struct containing the stepsizes for θ and ρ to use for the Metropolis-Hastings algorithm.
• nx::Int The number of x dimensions.
• ntheta::Int The number of θ dimensions.
• nloc::Int The number of unique settings of x in the experimental observations.

Details This function performs MCMC simulation to draw posterior samples for the parameters. First, the function determines whether to use a univariate data model or not. Then, samples are drawn from the posterior distributions using the following:

• Metropolis-Hastings updates for θ and ρ
• Gibbs updates for τ^2, σ^2, δ