Bayesian Uncertainty Quantification for Functional Response

Abstract

This chapter addresses the stochastic modeling of functional response, which is a major concern in engineering implementation. We first introduce a general framework and several conventional models for functional data, including the functional linear model, penalized regression splines, and the spatial temporal model. However, in engineering practice, a naive mathematical modeling of functional response may fail due to the lack of expressing the underlying physical mechanism. We propose a series of quasiphysical models to handle the functional response. A motivating example of metamaterial design is thoroughly discussed to demonstrate the idea of quasiphysical models. In real applications, various uncertainties have to be taken into account, such as that of the permittivity or permeability of the substrate of the metamaterial. For the propagation of uncertainty, simulation‐based methods are discussed. A Bayesian framework is presented to deal with the model calibration in the case of functional response. Experimental results illustrate the efficiency of the proposed method.

1. Introduction

In recent years, computer experiments have become widely adopted in both engineering applications and scientific research to replace or support their physical counterparts. Functional response is the mathematical representation of system behaviors, where the data are collected over an interval of some input indices. With the advance of modern simulation and experiment technology, accessing functional data becomes easier. Functional response can be in the form of one‐dimensional data such as a curve or higher dimensional data such as an image, which can provide better physical insights. However, even with the advancement of computer technology, full simulation based on a finite element method or a finite difference method still takes an extensive amount of time. To reduce the amount of simulation time, historical simulated data are usually used to build a cheaper metamodel [1], in which the functional response of unobserved input can be predicted by either regression or interpolation. The simplest representation of functional data can be considered basis expansion, where polynomials are used to formulate the input‐output relation [2]. For frequency response analysis, Fourier series are usually applied to replace the polynomials [3]. Both methods are categorized as linear regression, which requires parameter estimations. Nonparametric approaches were also used to analyze functional data in many scientific and engineering fields [4]. The purpose of building these models is to provide the “best” estimate regarding the given data, while providing a statistical scheme for prediction at unobserved inputs.

In this chapter, we provide a more sophisticated approach to naturally analyze functional responses, which may suggest more insightful conclusions that may not be apparent otherwise. We introduce one motivating example of functional response in computer experiments. In the design of metamaterial, the goal is to establish a relationship between the physical dimensions of a unit cell and its electromagnetic (EM) frequency response [5]. In practice, designers usually evaluate the EM properties of a metamaterial microstructure via full‐wave simulation data, such that corresponding adjustments are constantly made to the design (dielectric architecture, microstructure topologies, etc.) until a desired performance is achieved. Figure 1 depicts an example of unit cell design whose response phases differ on a frequency span along with the varying geometric parameter. Naïve regression‐based metamodels fail in dealing with such a problem because they require building regressions for each output, which could be very expensive and leaves the correlation between different frequencies unutilized. Moreover, when resonance is involved, the functional data cannot be well described by polynomials or splines. However, this can be overcome by some quasiphysical models, which explore the essential physical mechanism. In addition, a more general two‐stage modeling scheme can be applied, where in Stage I, we approximate the response with rational functions. This allows us to decompose the continuous response into a few discrete parameters. Stage II consists of a nonparametric metamodel to capture the input dependence.

Figure 1.

Example of functional response in metamaterial unit cell design‐phase shift between different physical dimension inputs.

2. General models for functional response

Various statistical models, including the spatial temporal model, functional linear model, and penalized regression splines, have been widely discussed in the past. Most models share a unified expression that sums up a mean function $\mu \text{(}f,x\text{)}$ and a random term $\epsilon \text{(}f,x\text{)}$, written as

$y\text{(}f,x\text{)}=\mu \text{(}f,x\text{)}+\epsilon \text{(}f,x\text{)}$E1

where $y$ is the response, $x=\{x_1,\mathrm{\dots },x_p\}$ is the input variables with dimensionality $p$, and $f$ represents some index, which could be the frequency of an electromagnetic wave or the time of a time series. Despite the shared form, these models differ in the way the mentioned terms are estimated.

2.1. Functional linear model

To model the functional response, the primary task is to estimate the mean function $\mu \text{(}f,x\text{)}$, on which a certain form is often imposed. As a generalization of linear regression models, the functional linear model is in the form of

$\mu \text{(}f,x\text{)}={\mathit{\beta }}_0\text{(}f\text{)}+x^T\mathit{\beta }\text{(}f\text{),}$E2

with basis functions $\mathit{\beta }={\displaystyle \left\{{\mathit{\beta }}_0\text{(}f\text{) ,}\beta \text{(}f\text{)}\right\}}$ ($\mathit{\beta }\text{(}f\text{)=(}{\mathit{\beta }}_1\text{(}f\text{),}\mathrm{\dots },{\mathit{\beta }}_p\text{(}f{\text{))}}^{{}^{\text{T}}}$ has the same dimensionality to the input variable), which incorporate the index dependence, and can be seen as an extension to the parameters of linear regression models. Therefore, by substituting the mean function into Eq. (1), we obtain the resulting output

$y\text{(}f,x\text{)}={\beta }_0\text{(}f\text{)}+x^T\mathit{\beta }\text{(}f\text{)}+\epsilon \text{(}f,x\text{)}\text{.}$E3

Given a certain index $f$, this model is a universal linear model. Furthermore, it contains an underlying index‐varying effect of x, whereas $\mathit{\beta }$ is assumed to be a smooth function of $f$. Thus, the model is referred to as a functional linear model. To estimate the coefficients ${\beta }_0\text{(}f\text{)}$ and $\beta \text{(}f\text{)}$, it is straightforward to apply the least squares method, which adopts the data collected at f. However, smoothing over f componentwise, using penalized splines, can enhance the efficiency of estimates [6].

Penalized regression splines implement estimation of smoothing basis functions in functional linear models by minimizing the penalized least squares. They are widely adopted in modeling functional responses due to their easy implementation and low computational cost [6].

Noted that the primary purpose of applying penalized regression splines is to estimate the basis function $\beta$. Suppose we have $n$ data ${\displaystyle \left\{{\displaystyle \left(f,y_i\right)},i=1,\mathrm{\dots },n\right\}}$, and the basis function is a random sample from

${\beta }_j=m\text{(}f\text{)}+{\delta }_j,$E4

with $j=1,\mathrm{\dots },p$. $m\text{(}f\text{)}$ is an unspecified smooth mean function of $\mathit{\beta }$ and ${\delta }_j$ is a zero mean random error. In practice, $m\text{(}f\text{)}$ can be estimated by a series of power‐truncated spline basis $1,f,f^2,f^p,{\text{(}f-{\kappa }_1\text{)}}_{+}^p,\dots ,{\text{(}f-{\kappa }_K\text{)}}_{+}^p$, where $\{{\kappa }_1,\dots ,{\kappa }_K\}$ is a given set of knots and $a_{+}$ denotes the positive part of $a$, i.e., $a_{+}=\text{(}a+\text{|}a\text{|})\text{/2}$. Therefore, the model in Eq. (4) can be approximately written as

${\beta }_j\approx {\alpha }_0+{\displaystyle \sum _{l=\text{1}}^p{\alpha }_l{f}^l}+{\displaystyle \sum _{k=\text{1}}^K{\alpha }_{k+p}{(f‐{\kappa }_k)}_{+}^p}+{\delta }_j,$E5

where ${\alpha }_j$ represents coefficients whose values can be obtained via least squares estimates. Generally, overfitting in the approximation of $m\text{(}f\text{)}$ may occur, which leads to high variance and poor prediction. To avoid large modeling bias, the trade‐off between model bias and overfitting requires careful consideration. In order to resolve such a problem, variable selection procedures should be applied to the linear regression model. However, when the number of involved basis functions is very large, variable selection would encounter great computational difficulty [7]. Alternatively, ${\alpha }_j$ is estimated by minimizing the penalized least squares function in the form of

${{\displaystyle \sum _{i=\text{1}}^n{\displaystyle \left[y_i-x_{ij}{\displaystyle \left\{{\alpha }_0+{\displaystyle \sum _{l=\text{1}}^l{\alpha }_l{f}^l}+{\displaystyle \sum _{k=\text{1}}^K{\alpha }_{k+p}{(f‐{\kappa }_k)}_{+}^p}\right\}}\right]}}}^2+g{\displaystyle \sum _{\text{k}=\text{1}}^K{\alpha }_{{}_{k+p}}^2},$E6

where $g$ is a tuning parameter determined by cross‐validation or generalized cross‐validation [8].

The smoothing method with penalized splines estimates also requires selection of the number of knots and the order p, which may vary from case to case. Fortunately, the estimates are not sensitive to these choices; and cubic splines are suggested in most cases [6], which ensure continuous second‐order and piecewise continuous third‐order derivatives at the knots. Meanwhile, knots are usually selected from the interval over which f is evenly distributed, or ${\kappa }_k$ is taken to be the $100k\text{/(}K+1\text{)th}$ percentile from the unevenly distributed f.

2.2. Spatial temporal model

The spatial temporal model is defined by the sum of a mean function, $\mu \text{(}f,x\text{)}$, and $\epsilon \text{(}f,x\text{)}$ of a zero‐mean Gaussian random field. It is a generalization of the Gaussian processes (GP) model, which has been widely adopted for spatial statistic problems [4, 9].

Both of the preceding models aim to represent the functional data in terms of their mean functions. In contrast, the spatial temporal model utilizes the property of the normal distribution of the residuals; thus, the output can be seen as a realization of a Gaussian random field. We assume a mean function $\mu \text{(}f,x\text{)}$ in the form of

$\mu \text{(}f,x\text{)}={\displaystyle \sum _{i=0}^n{h}_i\text{(}x\text{)}{\beta }_i\text{(}f\text{)}}\stackrel{def}{=}h{\text{(}x\text{)}}^T\mathit{\beta }\text{(}f\text{),}$E7

where $h\text{(}x\text{)}$ and $\mathit{\beta }\text{(}f\text{)}$ are two series of basis functions of the input variable and index variable, respectively. Such an assumption leads to a spatial temporal model

$y\text{(}f,x\text{)}=h{\text{(}x\text{)}}^T\mathit{\beta }\text{(}f\text{)}+\epsilon \text{(}f,x\text{),}$E8

where $\epsilon \text{(}f,x\text{)}$ is a zero‐mean Gaussian random field, and the covariance function follows the form

$cov\{\epsilon \text{(}f,x\text{)},\epsilon \text{(}f,x\text{'})\}=\text{K(}{\kappa }_f;\text{|}x‐x\text{'}\text{|),}$E9

where $K\text{(}{\kappa }_f\text{)}$ denotes the covariance matrix, whose (i,j) element $K\text{(}{\kappa }_f;|x_i-x_j|\text{)}$ measures the covariance between $x_i$ and $x_j$. ${\kappa }_f$ is an f‐dependent hyperparameter that controls the properties of the covariance.

Suppose we have obtained observation $y\text{(}{f}_j,x_i\text{)}$ at input sites $\text{(}{f}_j,x_i\text{)}$ with $j=1,\mathrm{\dots },J$ and $i=\text{1,}\mathrm{\dots },n$, where $J$ and $n$ are the length of indices and input settings. $\mathit{\beta }\text{(}f\text{)}$ and ${\kappa }_f$ can be calculated following the hyperparameter estimation procedure within a standard Gaussian processes model [4]. The spatial temporal model also allows predictions at unobserved sites $f_{*}$ and $x_{*}$. The procedures for prediction are summarized in the following algorithm.

Step 1: For $j=1,\mathrm{\dots },J$, calculate the best estimates of $\widehat{\mathit{\beta }}\text{(}f\text{)}$ and ${\kappa }_f$ for $f=f_j$ by maximizing the (log) likelihood, given by

$\text{log }p\text{(}y\text{|}X\text{)}=-\frac{1}{2}\text{[(}y-h\text{(}x\text{)}\mathit{\beta }{\text{)}}^T{K}^{-\text{1}}\text{(}{\kappa }_f\text{)(}y-h\text{(}x\text{)}\mathit{\beta }\text{)]}-\frac{1}{2}\text{log}\left|K\right|-\frac{N}{2}\text{log 2π}$E34

Step 2: According to the data $\{x_i,y\text{(}{f}_j,x_i)\}$, obtain estimates $\mu \text{(}f,x\text{)}$, and ${\widehat{\kappa }}_f$. Calculate prediction $y\text{(}{f}_j,x_{*}\text{)}$, at $f=f_j$, with the best linear unbiased prediction [2]

$y\text{(}{f}_j,x_{*}\text{)}=h{\text{(}{x}_{*}\text{)}}^T\widehat{\mathit{\beta }}+K^T\text{(}{\widehat{\kappa }}_f\text{;|}{x}_{*}‐x|\text{) }{K}^{‐1}\text{(}{\widehat{\kappa }}_f\text{;|}x-x\text{'}|\text{)(}y-h{\text{(}x\text{)}}^T\widehat{\mathit{\beta }}\text{)}$E35

Step 3: For new index $f_{*}\in \text{[}{f}_1,f_J\text{]}$ and given outputs at two existing indices $y\text{(}{f}_0\text{)}$ and $y\text{(}{f}_1\text{)}$, use linear interpolation to make predictions for y$\text{(}{f}_{*},x_{*}\text{)}$, as

$$y\text{(}{f}_{*},x_{*}\text{)}=y\text{(}{f}_0,x_{*}\text{)}+\text{(}{f}_{*}-f_0\text{)}\frac{y\text{(}{f}_1,x_{*}\text{)}-y\text{(}{f}_0,x_{*}\text{)}}{f_1-f_{*}}$$E30255

2.3. Quasiphysical model

Metamaterial frequency response, for example, modeling the resonance response is often quite challenging and cannot be achieved with the models introduced above. This is due to that the above models are based upon linear regression and simply encode the index dependence within the linear index‐dependent smooth basis functions. However, when distinct resonance peaks exist, a common scenario in radio frequency engineering, fitting to these smooth basis functions often, leads to poor accuracy [10]. To deal with these problems, we tend to utilize some underling physical mechanism and establish a quasiphysical modeling method. For example, the mean function $\mu \text{(}f\text{, }x\text{)}$ is represented by the combination of some link function $\mathbf{L(}f,•\text{)}$, which follows certain physical mechanisms, and a set of low‐dimensional scaling variables $\phi \text{(}x\text{)}$, i.e.,

$\mu {\displaystyle \left(f,x\right)}=\text{L}{\displaystyle \left(f,\phi \text{(}x\text{)}\right)}.$E10

Then, we have $y\text{(}f,x\text{)}=\text{L}{\displaystyle \left(f,\phi \text{(}x\text{)}\right)}+\epsilon \text{(}f,x\text{)}$. Instead of finding a single function with respect to both frequency index and input variables, the functional response is separated into two parts: a physical meaningful link function $\mathbf{L(}f,•\text{)}$ contains the functional features, whereas the other captures the relationship between input variables $x$ and scaling variables $\phi \text{(}x\text{)}$. This separation often leads to dimension reduction in statistical models. In the example of metamaterial design, the functional response is represented by the effective permittivity of a unit cell, which can be well fitted by a Drude‐Lorentz form [11],

$\mathbf{L(}f,\phi \text{(}x\text{)}\text{)}={\epsilon }_a{\displaystyle \left(1-\frac{F_e{f}^2}{f^2-f_0^2+i{\gamma }_{e}f}\right)}.$E11

where $\phi \text{(}x\text{)}\equiv \{{\epsilon }_a,F_e,f_0,{\gamma }_e\}$ is the intermediate variable which can be estimated via fitting the functional response by Eq. (11), meanwhile $\phi \text{(}x\text{)}$ is a function of input variables. Here, we choose the Gaussian processes (GP) regression model for interpolate new ${\phi }^{*}$ given previous obtained pairs $\{x,\phi \text{(}x\text{)}}$ and new $x^{*}$. Once the new ${\phi }^{*}$ is obtained, it can then be used to evaluate the new functional response by Eq. (11). Figure 2 displays a smooth surface of $f_{\text{0}}$ and an example of predicted effective permittivity.

Figure 2.

Example of modeling functional response assisted by the physical model: (a) Gaussian process surface of a scaling variable; (b) predicted functional response.

The aforementioned Drude‐Lorentz model allows high accuracy only when the metamaterial system works in a static or quasistatic regime, such that the metamaterial architecture can be seen as a single piece of effective medium. However, for complex metamaterial systems, the working regime is beyond static; thus, approximation accuracy by such a model is severely deteriorated. We noted that EM waves propagate through each layer of metamaterial like a current on transmission lines. Such a perspective transfers the EM field problems to circuit problems. Hence, function response to a continuous spectrum is reduced to discrete LRC (short form of inductor, resistor and capacitor) networks. We propose a two‐stage modeling scheme, where in the first stage, a vector fitting (VF) technique is adopted to provide accurate rational approximation to frequency responses with distinct resonances. Its results are easily interpreted as an equivalent circuit. The approximation accuracy to a frequency response and its corresponding equivalent circuits are shown in part (a) and (b) of Figure 3, respectively. And in Stage II, the empirical circuit elements are then taken as the target response in statistical models to establish the mapping input‐output relation by performing regression, which also allows predictions at unobserved input sites. Part (c) of Figure 3 presents the GP surface built of circuit parameters over two input variables. A graphical display of this two‐stage approach is illustrated in Figure 4. To predict functional response at unobserved input, it is implemented by first predict the presenting circuit elements and then recover the response.

Figure 3.

Example of modeling frequency response via (a) vector fitting, (b) equivalent circuit, and (c) GP regression.

Figure 4.

Flowchart of the two‐stage modeling approach.

3. Uncertainty quantification

In the engineering modeling and design, uncertainty is ubiquitous, due to the inability to specify a “true” input or model parameter. Quantifying the uncertainty of the model, e.g., in the form of predictive confidence intervals, is of great importance for decision making and advanced design [1]. In general, uncertainty quantification can be divided into two major types of problem: forward uncertainty propagation and inverse assessment of model and parameter uncertainties [12].

The full relationship between experimental output $z\text{(}f,x\text{)}$ and simulation output $y\text{(}f,x\text{)}$ can be expressed as

$z\text{(}f,x)=y\text{(}f,x\text{)}+{\epsilon }^{\prime }(f)=\text{L}{\displaystyle \left(f,\phi \text{(}x,\theta \text{)}+\eta \text{(}x\text{)}\right)}+\epsilon \text{(}f,x\text{)+}{\epsilon }^{\prime }\text{(}f\text{),}$E12

where $\phi \text{(}x,\mathit{\theta }\text{)}$ denotes the GP regression model, which depends on the input variable $x$ and several unobservable calibration parameters $\theta$. $\eta {\displaystyle \left(x\right)}$ is the additive discrepancy function (or model bias function), which does not depend upon the calibration parameters. To reduce the complexity of the analysis, we assume that $\epsilon \text{(}f,x\text{)}$ is a zero‐mean Gaussian random field and being independent of $x$. And then we can merge $\epsilon \text{(}f,x\text{)}$ and $\epsilon \text{'}\text{(}f\text{)}$ together, which is denoted by $\epsilon (f)$. Thus, the model becomes

$z\text{(}f,x)=\text{L}{\displaystyle \left(f,\phi \text{(}x,\theta \text{)}+\eta \text{(}x\text{)}\right)}+\epsilon \text{(}f\text{)}$E13

where $\epsilon (f)$ is assumed to be a zero‐mean Gaussian noise with known variance $\lambda$, $\epsilon (f_i)\sim \mathcal{N}(0,\lambda )$.

In forward problems, with an uncertain input x and given model parameters θ, the model output $y$ and other quantities of interest are to be calculated. On the other hand, the inverse problem is to estimate the values of model parameters θ such that it makes the model’s output fit the experimental data as accurate as possible (or satisfy some precision requirements).

3.1. Metamodel‐based uncertainty propagation

The main problem in analyzing uncertainty propagation is obtaining an analytical representation of the metamodel for any arbitrary (uncertain) input values. Given its probability density, the Bayesian framework can provide a probability measure of random inputs on the output field. The purpose of such an operation is to evaluate the influence of an uncertain input on the model response.

Assume that the Gaussian process regression model is trained on a dataset with the input $X={\displaystyle \left\{x_1,\mathrm{\dots },x_N\right\}}$ and the corresponding intermediate variable $\Psi ={{\displaystyle \left(\phi \text{(}{x}_1\text{)},\mathrm{\dots },\phi \text{(}{x}_N\text{)}\right)}}^{\text{T}}$ which is obtained by fitting algorithm in Stage I. The GP hyperparameters learned from the data are denoted by $\gamma$. The uncertainty of the input variable $x^{*}$ is captured by a probability density function,

$x^{*}\sim p\text{(}{x}^{*}\text{)}$E14

At a deterministic test input $x^{\text{*}}$, the predictive distribution of the function, $p\text{(}{\phi }^{\text{*}}{\text{|x}}^{\text{*}},X,\Psi ,\gamma \text{)}$ (for simplicity, we use ${\phi }^{\text{*}}$ to denote $\phi \text{(}{x}^{\text{*}}\text{)}$, the output of the metamodel.), is Gaussian with mean

${\tilde{\phi }}^{\text{*}}=\mathbf{E}\text{(}{\phi }^{\text{*}}{\text{|x}}^{\text{*}},X,\Psi ,\gamma \text{)}={\displaystyle \sum _{i=1}^N{\zeta }_i}\text{C(}{x}_i,{\text{x}}^{\text{*}}\text{),}$E15

and variance

$cov\text{(}{\varphi }^{\text{*}}\text{)}=C{\text{(x}}^{\text{*}},{\text{x}}^{\text{*}}\text{)}-{\displaystyle \sum _{i,j=1}^N\text{(}C-{\sigma }^{2}I}{\text{)}}^{-1}{\text{C(x}}^{\text{*}},x_i\text{)}C{\text{(x}}^{\text{*}},x_j\text{)}$E16

where ${\zeta }_i$ is the ith element of column vector $\zeta =\text{[}C\text{+}{\sigma }^{2}I\text{]}\Psi$. $C$ denotes the covariance matrix of the Gaussian process, whose ijth element is given by $C_{ij}=C\text{(}{x}_i,x_j\text{)}$.

The final goal is to propagate uncertainty through the link function $\text{L}(f,{\phi }^{*})$. The computation of the statistics is implemented by integrating over the uncertainty with the mean

${\mu }_{{\text{L}}^{\text{*}}}={\displaystyle \int \text{L(}f,{\phi }^{*}\text{)}p\text{(}{\phi }^{*}{\text{|x}}^{\text{*}},X,\Psi ,\gamma \text{)}p}\text{(}{x}^{*}\text{)}d{\phi }^{*}d{x}^{*}.$E17

and the variance

${\sigma }_{{\text{L}}^{*}}^2={\displaystyle \int {\text{[L(}f,{\phi }^{*}-{\mu }_{{\text{L}}^{*}}\text{)]}}^{2}p\text{(}{\phi }^{*}\text{|}{x}^{\text{*}},X,\Psi ,\gamma \text{)}p}\text{(}{x}^{*}\text{)}d{\phi }^{*}d{x}^{*}.$E18

The uncertainty propagation is induced by the variability of the input variable. For example, in metamaterial engineering, the dimension of a design parameter, say the thickness of the metallic microstructure layer, could differ from what has been instructed during the manufacturing processes. From measurements, the value of such a variable would rather follow a distribution than be pre‐specified as an exact value. Therefore, the analysis of uncertainty propagation is needed to be in the metamaterial design process.

3.2. Bayesian calibration

Compared to uncertainty forward propagation, the inverse problem is more difficult yet of great importance in enhancing the fidelity of metamodels. Two major aspects concerning the inverse problem are measuring model discrepancy and model calibration. In this chapter, we use the formulation to address both issues within an updating process, similar to that proposed in Ref. [12].

3.2.1. The model

In this section, we introduce the details of performing Bayesian calibration with regard to Eq. (13). The calibration parameters, denoted by $\mathit{\theta }$, are defined as any physical parameters that can be specified as an input to the statistical model given by Eq. (13). The fundamental difference between $x$ and $\mathit{\theta }$ is that the former refers to design inputs whose value can be specified by the user during experiment and simulation, whereas the latter cannot be controlled and its true value is not directly observable [12]. In the previous chapter, the calibration parameter is not explicitly specified. However, we here include it in the framework to quantify its uncertainty, which completes the full cycle of metamaterial design and modeling. Suppose $\mathit{\theta }$ represent a constitutive parameter, say permittivity, of a dielectric used to fabricate the metamaterial system, which cannot be accurately measured directly.

Before offering the detailed statistics for uncertainty quantification, we must note that the purpose of parameter calibration is to provide an accurate prediction with the metamodel with a small amount of data. An even smaller amount of experimental data is acquired to calibrate and validate the main model. To select the “best” experiment samples, uniform experimental design techniques are usually applied [6]. A Latin hypercube sampling, for example, is widely used for such cases, mainly due to its good coverage property [13].

The data corresponding to the metamodel $\phi$ are obtained at $D_1={\displaystyle \left\{\text{(}{x}_1^{\text{'}},{\theta }_1\text{),}\mathrm{\dots }\text{,(}{x}_N^{\text{'}},{\theta }_N\text{)}\right\}}$, where $\{x_1^{\text{'}},\mathrm{\dots },x_N^{\text{'}}\}$ and $\{{\theta }_1,\mathrm{\dots },{\theta }_N\}$ are the set of design inputs and calibration parameters. Although the notation is included, the true values of the calibration parameters are unknown throughout the entire calibration process. The inverse problem of uncertainty quantification is implemented in an updated formulation with a Bayesian approach [1]. In model (13), the metamodel, $\mathit{\phi }\text{(}x,\mathit{\theta }\text{)}$, and discrepancy function, $\eta \text{(}x\text{)}$, are both Gaussian processes:

$\phi \text{(}x,\theta )\sim \mathcal{N}{\displaystyle \left(m_1\text{(}x,\theta \text{),}{C}_1{\displaystyle \left(\text{(}x,\theta \text{)},\text{(}x\text{'},\theta \text{'}\text{)}\right)}\right)},$E19

$\eta \text{(}x)\sim \mathcal{N}{\displaystyle \left(m_2\text{(}x\text{),}{C}_2\text{(}x,x\text{'}\text{)}\right)},$E20

where $m_1\text{(}x,\theta )=h_1{\text{(}x,\theta \text{)}}^{\text{T}}{\beta }_1$ and $m_2\text{(}x)=h_2{\text{(}x\text{)}}^{\text{T}}{\beta }_2$ [12]. $C_1{\displaystyle \left({\displaystyle \left(\cdot ,\cdot \right)},{\displaystyle \left(\cdot ,\cdot \right)}\right)}$ and $C_2{\displaystyle \left(\cdot ,\cdot \right)}$ are covariance functions, which can be parameterized by some hyperparameters, denoted by ${\Gamma }_{\text{1}}$ and ${\Gamma }_{\text{2}}$, respectively. Let us denote these hyperparameters with $\Gamma =\text{(}{\Gamma }_{\text{1}},{\Gamma }_{\text{2}}\text{)}$, collectively. There are many candidates of covariance functions from which one can chose. For example, as one of the most applied covariance functions, squared exponential function, in form of

$\begin{array}{c}{C}_1{\displaystyle \left({\displaystyle \left(x,\theta \right)},{\displaystyle \left(x\text{'},\theta \text{'}\right)}\right)}={\sigma }_1^2\text{exp}{\displaystyle \left\{-{\text{(}x-x\text{'}\text{)}}^{\text{T}}{V}_{\text{1x}}\text{(}x-x\text{'}\text{)}\right\}}\text{exp}{\displaystyle \left\{-{\text{(}\theta -\theta \text{'}\text{)}}^{\text{T}}{V}_{\theta }\text{(}\theta -\theta \text{'}\text{)}\right\}},\\ C_2{\displaystyle \left(x,x\text{'}\right)}={\sigma }_2^2\text{exp}{\displaystyle \left\{-{\text{(}x-x\text{'}\text{)}}^{\text{T}}{V}_{\text{2x}}\text{(}x-x\text{'}\text{)}\right\}}\text{ ,}\end{array}$E21

can provide smooth samples to infer the latent function variable. In Eq. (21), the value of hyperparameters can be inferred via Markov Chain Monte Carlo (MCMC) techniques.

3.2.2. Data and prior distribution

Let us denote the matrix of basis functions $H_1\text{(}{D}_{\text{1}}\text{)}$ with rows ${\displaystyle \left\{h_1{\text{(}{x}_1^{\text{'}},{\theta }_1\text{)}}^T,\mathrm{\dots },h_1{\text{(}{x}_N^{\text{'}},{\theta }_N\text{)}}^T\right\}}$, which leads to the expectation of $\phi$ as $H_1\text{(}{D}_{\text{1}}\text{)}{\beta }_1$. Similarly, from the experimental observations we can obtain $\widehat{\phi }$, the estimation of $\phi$ by Eq. (13). It can be further augmented by the calibration parameter at each $x$, with $D_2\text{(}\theta \text{)}={\displaystyle \left\{\text{(}{x}_1,\mathit{\theta }\text{)},\mathrm{\dots },\text{(}{x}_n,\mathit{\theta }\text{)}\right\}}$. In contrast to the simulation output $\{x_1^{\text{'}},\mathrm{\dots },x_N^{\text{'}}\}$, the experimental data are usually acquired with much smaller size, i.e., $n<<N$, which is in accordance with the purpose of reducing the amount of physical experiments with calibrated models. Meanwhile, we use $x_i$ and $x_i\text{'}$ to describe that the observation points could be different between two datasets. The expectation $\widehat{\phi }$ of can be represented by $H_1{\displaystyle \left\{D_{\text{2}}\text{(}\theta \text{)}\right\}}{\beta }_1+H_2\text{(}{D}_{\text{2}}\text{)}{\beta }_2$. We write the full data vector ${\Omega }^T=\{{\phi }^T,{\widehat{\phi }}^T\}$, which is obtained via Stage I given the simulation and observation of functional response. Meanwhile, they are normally distributed given the full set of parameters $\{\mathit{\theta },\mathit{\beta },\mathit{\varphi }\}$ ($\beta ={\text{(}{\beta }_1^T,{\beta }_2^T\text{)}}^T$, $\mathit{\varphi }=(\lambda ,\Gamma \text{)}$).

The goal of calibration is to obtain $p\text{(}\mathit{\theta }|\Omega \text{)}$, the posterior distribution of conditional only on the full data $\Omega$. To derive the posterior distribution of parameters, we begin with the normal distribution of the full set of data, during which the likelihood function will yield a Gaussian [14], with mean

$E\text{(}\Omega |\mathit{\theta },\mathit{\beta },\mathit{\phi }\text{)}=m_d\text{(}\mathit{\theta })=H\text{(}\mathit{\theta }\text{)}\mathit{\beta },$E22

where

$H\text{(}\theta \text{)}={\displaystyle \left(\begin{array}{cc}{H}_1\text{(}{D}_{\text{1}}\text{)}& 0\\ H_1{\displaystyle \left\{D_{\text{2}}\text{(}\theta \text{)}\right\}}& H_2\text{(}{D}_{\text{2}}\text{)}\end{array}\right)}.$E23

To specify the variance matrix of $\Omega$, we need the variance matrix of $\phi$, denoted by $V_1\text{(}{D}_{\text{1}}\text{)}$, whose (i,i') element is $C_1{\displaystyle \left(\text{(}{x}_i^{\text{'}}\text{ ,}{\theta }_i\text{),(}{x}_{i\text{'}}^{\text{'}},{\theta }_{i\text{'}}\text{)}\right)}$. Similarly, we can define $V_1{\displaystyle \left\{D_2\text{(}\mathit{\theta }\text{)}\right\}}$ and $V_2\text{(}{D}_2\text{)}$. Let $C_1{\displaystyle \left\{D_1,D_2\text{(}\mathit{\theta }\text{)}\right\}}$ be the matrix with (i,j) element $C_1{\displaystyle \left\{\text{(}{x}_i^{\text{'}}\text{ ,}{\theta }_i\text{),(}{x}_j^{\text{'}},{\theta }_j\text{)}\right\}}$. Therefore,

$Var\text{(}\Omega \text{|}\mathit{\theta },\mathit{\beta },\mathit{\phi }\text{)}=V_d\text{(}\mathit{\theta }\text{)}={\displaystyle \left(\begin{array}{cc}{V}_1\text{(}{D}_{\text{1}}\text{)}& C_1{{\displaystyle \left\{D_1,D_2\text{(}\mathit{\theta }\text{)}\right\}}}^T\\ C_1{\displaystyle \left\{D_1,D_2\text{(}\mathit{\theta }\text{)}\right\}}& \lambda I\text{+}{V}_1{\displaystyle \left\{D_2\text{(}\mathit{\theta }\text{)}\right\}}\text{+}{V}_2\text{(}{D}_2\text{)}\end{array}\right)},$E24

where $I$ is the $n\times n$ identity matrix.

To derive the posterior distribution under the Bayesian framework, the prior distributions of parameters, $\{\mathit{\theta },\mathit{\beta },\mathit{\varphi }\}$, must also be independently specified. Following the suggestion of [12], we chose conjugate prior for $\mathit{\theta }$ and $\mathit{\varphi }$, and a weak prior for $\beta$, specifically $p\text{(}{\beta }_1,{\beta }_2\text{)}\propto \text{1}$, then we have $p\text{(}\mathit{\theta },\mathit{\beta },\mathit{\varphi }\text{)}=p\text{(}\mathit{\theta }\text{)}p\text{(}\mathit{\varphi }\text{)}$, where $p\text{(}\mathit{\varphi })=p\text{(}\lambda \text{)}p\text{(}{\Gamma }_1\text{)}p\text{(}{\Gamma }_2\text{)}$. Meanwhile, Bayesian inference with MCMC requires specification of proper prior distributions to perform Bayesian statistics. For such purpose, conjugate priors are specified, e.g.

$\begin{array}{l}{\sigma }_1^2,{\sigma }_2^2\sim \text{I}\mathcal{G}\text{(a,b),}\\ V_{\text{1x}},V_{\text{2x}}\sim \mathcal{W}\text{(}\rho ,\nu \text{),}\\ \theta \sim \mathcal{N}\text{(}{\mu }_{\theta }{\text{,V}}_{\theta }\text{)}\text{.}\end{array}$E25

where $\text{I}\mathcal{G}$, $\mathcal{W}$, and $\mathcal{N}$ are inverse gamma, Wishart, and normal distributions, respectively [15].

3.2.3. Posterior distribution

Conditional on full data, the independence of parameters leads to the full joint posterior distribution

$\begin{array}{c}p\text{(}\mathit{\theta },\mathit{\beta },\mathit{\varphi }|\Omega \text{)}\propto p\text{(}\mathit{\theta }\text{)}p\text{(}\mathit{\varphi }\text{)}p\text{(}\mathit{\beta }\text{)}p\text{(}\Omega |m_d\text{(}\mathit{\theta }\text{),}{V}_d\text{(}\mathit{\theta }\text{))}\\ \propto p\text{(}\mathit{\theta }\text{)}p\text{(}\mathit{\varphi }\text{)}{\left|V_d^{-1}\text{(}\mathit{\theta }\text{)}\right|}^{-1/2}\\ \times \text{exp}{\displaystyle \left[-\frac{1}{2}{\displaystyle \left\{{{\displaystyle \left(\Omega -m_d\text{(}\mathit{\theta }\text{)}\right)}}^T{V}_d^{-1}\text{(}\mathit{\theta }\text{)}{\displaystyle \left(\Omega -m_d\text{(}\mathit{\theta }\text{)}\right)}\right\}}\right]}.\end{array}$E26

To obtain $p\text{(}\mathit{\theta }|\Omega \text{)}$, it is required to integrate out $\beta$ and hyperparameters $\mathit{\varphi }$ from Eq. (26). Integrating $\beta$ yields

$\begin{array}{c}p\text{(}\mathit{\theta },\mathit{\varphi }|\Omega \text{)}\propto p\text{(}\mathit{\theta }\text{)}p\text{(}\mathit{\varphi }\text{)}{\left|V_d^{-1}\text{(}\mathit{\theta }\text{)}\right|}^{-1/2}{\left|W\text{(}\mathit{\theta }\text{)}\right|}^{1/2}\hfill \\ \times \text{exp}{\displaystyle \left[-\frac{1}{2}{\displaystyle \left\{{{\displaystyle \left(\Omega -H\text{(}\mathit{\theta }\text{)}\widehat{\beta }\text{(}\mathit{\theta }\text{)}\right)}}^T{V}_d^{-1}\text{(}\mathit{\theta }\text{)}{\displaystyle \left(\Omega -H\text{(}\mathit{\theta }\text{)}\widehat{\beta }\text{(}\mathit{\theta }\text{)}\right)}\right\}}\right]},\hfill \end{array}$E27

where

$\widehat{\mathit{\beta }}\text{(}\mathit{\theta }\text{)}=W\text{(}\mathit{\theta }\text{)}H{\text{(}\mathit{\theta }\text{)}}^T{V}_d{\text{(}\mathit{\theta }\text{)}}^{‐1}\Omega ,$E28

$W\text{(}\mathit{\theta }\text{)}={{\displaystyle \left(H{\text{(}\mathit{\theta }\text{)}}^T{V}_d{\text{(}\mathit{\theta }\text{)}}^{‐1}H\text{(}\mathit{\theta }\text{)}\right)}}^{-1}.$E29

3.2.4. Calibration and prediction

Since the posterior distribution specified in Eq. (27) is a highly intractable function of $\mathit{\varphi }$, we need Monte Carlo method to integrate out $\mathit{\varphi }$ and get the numerical estimation for the posterior distribution of the calibration parameters $\theta$. The formulation is given by

$p\text{(}\theta \text{|}\Omega \text{)}=\frac{1}{M}{\displaystyle \sum _{i=1}^{M}p\text{(}\theta ,{\mathit{\varphi }}^i\text{|}\Omega \text{)}}.$E30

However, the purpose of calibration of parameters is to predict the real process rather than achieve their values. Therefore, in practice, we are rather more interested in expressing the posterior distribution of $\mathit{\varphi }$, which is a Gaussian process as well, conditional on the calibration parameters and estimated hyperparameters. The mean and covariance function of this GP is given by

$E{\displaystyle \left(\mathit{\phi }\text{(}x\text{)|}\mathit{\theta },\mathit{\varphi },\Omega \right)}=h{\text{(}x,\mathit{\theta }\text{)}}^T\widehat{\mathit{\beta }}\text{(}\mathit{\theta }\text{)+}t{\text{(}x,\mathit{\theta }\text{)}}^T{V}_d^{-1}\text{(}\mathit{\theta }\text{)}{\displaystyle \left(\Omega -H\text{(}\mathit{\theta }\text{)}\widehat{\mathit{\beta }}\text{(}\mathit{\theta }\text{)}\right)},$E31

where $h\text{(}x,\theta \text{)}={\displaystyle \left(\begin{array}{c}{h}_1\text{(}x,\mathit{\theta })\\ h_2\text{(}x\text{)}\end{array}\right)},t\text{(}x,\mathit{\theta })={\displaystyle \left(\begin{array}{c}{V}_1{\displaystyle \left(\text{(}x,\mathit{\theta }\text{),}{D}_1\right)}\\ V_1{\displaystyle \left(\text{(}x,\mathit{\theta }\text{),}{D}_2\text{(}\mathit{\theta }\text{)}\right)}+V_2{\displaystyle \left(x,D_2\right)}\end{array}\right)}$,

and covariance

$\begin{array}{l}cov{\displaystyle \left(\phi \text{(}x\text{),}\phi \text{(}x\text{'}\text{)|}\mathit{\theta },\mathit{\varphi },\Omega \right)}=c_1{\displaystyle \left(\text{(}x,\mathit{\theta }\text{)},\text{(}x\text{'},\mathit{\theta }\text{)}\right)}+c_2\text{(}x,x\text{'}\text{)}-t\text{(}x,\theta {\text{)}}^T{V}_d^{-1}\text{(}\mathit{\theta }\text{)}t\text{(}x\text{'},\mathit{\theta }\text{)}\\ \text{+}{{\displaystyle \left(h\text{(}x,\mathit{\theta }\text{)-}H{\text{(}\mathit{\theta }\text{)}}^T{V}_d^{-1}\text{(}\mathit{\theta }\text{)}t\text{(}x,\mathit{\theta }\text{)}\right)}}^{T}W\text{(}\mathit{\theta }\text{)}{\displaystyle \left(h\text{(}x\text{'},\mathit{\theta }\text{)-}H{\text{(}\mathit{\theta }\text{)}}^T{V}_d^{-1}\text{(}\mathit{\theta }\text{)}t\text{(}x\text{'},\mathit{\theta }\text{)}\right)}.\end{array}$E32

Inference about $\mathit{\varphi }\text{(}x\text{)}$ can be implemented again numerically with its posterior mean $E{\displaystyle \left[\phi \text{(}x\text{)|}\mathit{\theta },\mathit{\varphi },\Omega \right]}$ at estimated $\mathit{\theta }$ and $\mathit{\varphi }$, by integrating Eq. (31) with regard to Eq. (28). Given the estimation of $\varphi \text{(}x\text{)}$, the analysis of $\text{z}$ becomes straightforward by applying the link function $\text{L}{\displaystyle (\cdot )}$ as described in model (13).

So far, we have accomplished calibrating a metamodel in the Bayesian framework using the experimental data, which accounts for parameter uncertainty and corrects the model discrepancy and experimental uncertainty.

4. Simulation study

This section demonstrate the results obtained using the Bayesian uncertainty quantification framework for the metamaterial design problem with the models described in Sections 2 and 3, with examples. Of both propagation and inverse assessment, the overall model is formulated in Eq. (13), where geometric variable $w$ and incident angle $\alpha$ are input variables specified in the simulation, i.e., $x^T={\{w,\alpha \}}^T$. Thus, the model is expressed as

$\begin{array}{c}z\text{(}x)=y\text{(}x,\mathit{\theta }\text{)+}\eta \text{(}x\text{)+}\epsilon \hfill \\ =L{\displaystyle \left(f,\phi {\displaystyle \left(\{w,\alpha \},\mathit{\theta }\right)}+\eta {\displaystyle \left(w,\alpha \right)}\right)}+{\epsilon }^{\prime },\end{array}$E33

To demonstrate parameter calibration within the metamaterial modeling and design, we consider an example where the real part of the permittivity of a dielectric material, ${\epsilon }_d$, is defined as the calibration parameters $\theta ={\epsilon }_d$, and its prior is given normal distribution as model (25), with mean ${\mu }_{\theta }\text{=3}$ and variance ${\text{V}}_{\theta }\text{=0}\text{.5}$. Figure 5 illustrates the probability density function of this prior distribution. We demonstrate a measure of uncertain propagation in Figure 6, where predictions with 2.5% quantile (green or light gray) and 97.5% quantile (red or dark gray) of the samples are depicted to show the discrepancy induced by the uncertain input. Following the methodology introduced in Section 3, metamodels can be established for the simulation data and discrepancy function, with Gaussian process regression models. In our example, we obtained 92 simulation data to build GPs and 20 observations for calibration. The posterior distribution of the calibration parameter is also displayed in Figure 5. After calibration, the distribution of calibration parameter has a much smaller variance. The comparison between the prediction at posterior mean (cyan curve) and “real data” (blue dash) is shown in Figure 6, where the discrepancy reduction is remarkable.

Figure 5.

Comparison of prior and posterior distributions of the calibration parameter. The mean of the Gaussian distribution shifts from 3 to 3.17, and the variance is much smaller after Bayesian calibration.

Figure 6.

The effect of uncertainty propagation and results of parameter calibration.

5. Conclusion

In this chapter, we review several conventional model for functional response and present the quasiphysical model for functional response. Compared with the conventional models, this model can reveal the physical insight more clearly and make better use of historical experience. The two‐stage method was presented to model the frequency response of metamaterial and facilitate the design process. Using this approach, we decomposed the complex modeling problem into a vector fitting‐based equivalent circuit modeling process and a GP regression process, which can easily generate the mapping function from the structure’s geometric design. The predictive property of this model enables the massive reduction of time‐consuming simulations.

Another important topic with this chapter was the development and application of a Bayesian uncertainty quantification approach in dealing with functional response. Both forward uncertainty propagation and inverse assessment of the model were discussed, and a Bayesian framework was presented with simulation experimental results to deal with the model calibration for functional response. We envision that our two‐stage approach can be generalized to model any functional responses of a rational form. With the Bayesian framework for the functional data of computer experiments, we were able to incorporate our prior knowledge into the model and obtain a probabilistic measure of the uncertainty associated with metamaterial system design. This general methodology enables researchers and designers to achieve high efficiency and accuracy in modeling functional response with a considerably small amount of data. With a Bayesian calibration framework, we are able to constantly increase the precision of predictions of the functional response at unobserved sites, thus replacing expensive physical experiments.

Acknowledgments

This work was supported by the Shenzhen Key Laboratory of Data Science and Modeling.