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.
- functional response
- meta model
- Bayesian uncertainty quantification
- model calibration
- metamaterial design
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 , 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 . For frequency response analysis, Fourier series are usually applied to replace the polynomials . 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 . 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 . 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.
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 and a random term , written as
where is the response, is the input variables with dimensionality , and 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 , on which a certain form is often imposed. As a generalization of linear regression models, the functional linear model is in the form of
with basis functions (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
Given a certain index , this model is a universal linear model. Furthermore, it contains an underlying index‐varying effect of x, whereas is assumed to be a smooth function of . Thus, the model is referred to as a functional linear model. To estimate the coefficients and , 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 .
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 .
Noted that the primary purpose of applying penalized regression splines is to estimate the basis function . Suppose we have data , and the basis function is a random sample from
with . is an unspecified smooth mean function of and is a zero mean random error. In practice, can be estimated by a series of power‐truncated spline basis , where is a given set of knots and denotes the positive part of , i.e., . Therefore, the model in Eq. (4) can be approximately written as
where represents coefficients whose values can be obtained via least squares estimates. Generally, overfitting in the approximation of 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 . Alternatively, is estimated by minimizing the penalized least squares function in the form of
where is a tuning parameter determined by cross‐validation or generalized cross‐validation .
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 , 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 is taken to be the percentile from the unevenly distributed f.
2.2. Spatial temporal model
The spatial temporal model is defined by the sum of a mean function, , and 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 in the form of
where and are two series of basis functions of the input variable and index variable, respectively. Such an assumption leads to a spatial temporal model
where is a zero‐mean Gaussian random field, and the covariance function follows the form
where denotes the covariance matrix, whose (i,j) element measures the covariance between and . is an f‐dependent hyperparameter that controls the properties of the covariance.
Suppose we have obtained observation at input sites with and , where and are the length of indices and input settings. and can be calculated following the hyperparameter estimation procedure within a standard Gaussian processes model . The spatial temporal model also allows predictions at unobserved sites and . The procedures for prediction are summarized in the following algorithm.
Step 1: For , calculate the best estimates of and for by maximizing the (log) likelihood, given by
Step 2: According to the data , obtain estimates , and . Calculate prediction , at , with the best linear unbiased prediction 
Step 3: For new index and given outputs at two existing indices and , use linear interpolation to make predictions for y, as
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 . To deal with these problems, we tend to utilize some underling physical mechanism and establish a quasiphysical modeling method. For example, the mean function is represented by the combination of some link function , which follows certain physical mechanisms, and a set of low‐dimensional scaling variables , i.e.,
Then, we have . 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 contains the functional features, whereas the other captures the relationship between input variables and scaling variables . 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 ,
where is the intermediate variable which can be estimated via fitting the functional response by Eq. (11), meanwhile is a function of input variables. Here, we choose the Gaussian processes (GP) regression model for interpolate new given previous obtained pairs and new . Once the new is obtained, it can then be used to evaluate the new functional response by Eq. (11). Figure 2 displays a smooth surface of and an example of predicted effective permittivity.
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.
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 . In general, uncertainty quantification can be divided into two major types of problem: forward uncertainty propagation and inverse assessment of model and parameter uncertainties .
The full relationship between experimental output and simulation output can be expressed as
where denotes the GP regression model, which depends on the input variable and several unobservable calibration parameters . 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 is a zero‐mean Gaussian random field and being independent of . And then we can merge and together, which is denoted by . Thus, the model becomes
where is assumed to be a zero‐mean Gaussian noise with known variance , .
In forward problems, with an uncertain input x and given model parameters θ, the model output 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 and the corresponding intermediate variable which is obtained by fitting algorithm in Stage I. The GP hyperparameters learned from the data are denoted by . The uncertainty of the input variable is captured by a probability density function,
At a deterministic test input , the predictive distribution of the function, (for simplicity, we use to denote , the output of the metamodel.), is Gaussian with mean
where is the ith element of column vector . denotes the covariance matrix of the Gaussian process, whose ijth element is given by .
The final goal is to propagate uncertainty through the link function . The computation of the statistics is implemented by integrating over the uncertainty with the mean
and the variance
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. .
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 , 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 and 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 . 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 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 . A Latin hypercube sampling, for example, is widely used for such cases, mainly due to its good coverage property .
The data corresponding to the metamodel are obtained at , where and 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 . In model (13), the metamodel, , and discrepancy function, , are both Gaussian processes:
where and . and are covariance functions, which can be parameterized by some hyperparameters, denoted by and , respectively. Let us denote these hyperparameters with , 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
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 with rows , which leads to the expectation of as . Similarly, from the experimental observations we can obtain , the estimation of by Eq. (13). It can be further augmented by the calibration parameter at each , with . In contrast to the simulation output , the experimental data are usually acquired with much smaller size, i.e., , which is in accordance with the purpose of reducing the amount of physical experiments with calibrated models. Meanwhile, we use and to describe that the observation points could be different between two datasets. The expectation of can be represented by . We write the full data vector , 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 (, ).
The goal of calibration is to obtain , the posterior distribution of conditional only on the full data . 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 , with mean
To specify the variance matrix of , we need the variance matrix of , denoted by , whose (i,i') element is . Similarly, we can define and . Let be the matrix with (i,j) element . Therefore,
where is the identity matrix.
To derive the posterior distribution under the Bayesian framework, the prior distributions of parameters, , must also be independently specified. Following the suggestion of , we chose conjugate prior for and , and a weak prior for , specifically , then we have , where . Meanwhile, Bayesian inference with MCMC requires specification of proper prior distributions to perform Bayesian statistics. For such purpose, conjugate priors are specified, e.g.
where , , and are inverse gamma, Wishart, and normal distributions, respectively .
3.2.3. Posterior distribution
Conditional on full data, the independence of parameters leads to the full joint posterior distribution
To obtain , it is required to integrate out and hyperparameters from Eq. (26). Integrating yields
3.2.4. Calibration and prediction
Since the posterior distribution specified in Eq. (27) is a highly intractable function of , we need Monte Carlo method to integrate out and get the numerical estimation for the posterior distribution of the calibration parameters . The formulation is given by
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 , 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
Inference about can be implemented again numerically with its posterior mean at estimated and , by integrating Eq. (31) with regard to Eq. (28). Given the estimation of , the analysis of becomes straightforward by applying the link function 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 and incident angle are input variables specified in the simulation, i.e., . Thus, the model is expressed as
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, , is defined as the calibration parameters , and its prior is given normal distribution as model (25), with mean and variance . 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.
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.
This work was supported by the Shenzhen Key Laboratory of Data Science and Modeling.