{'fraction_non_alphanumeric': 0.0392046369344711, 'fraction_numerical': 0.010243922073740138, 'mean_word_length': 4.064621343390073, 'pattern_counts': {'":': 0, '<': 0, '<?xml version=': 0, '>': 1, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Hierarchical Stochastic Model in Bayesian Inference: Theoretical Implications and Efficient Approximation

We classify two types of Hierarchical Bayesian Model found in the literature as Hierarchical Prior Model (HPM) and Hierarchical Stochastic Model (HSM). Then, we focus on studying the theoretical implications of the HSM. Using examples of polynomial functions, we show that the HSM is capable of separating different types of uncertainties in a system and quantifying uncertainty of reduced order models under the Bayesian model class selection framework. To tackle the huge computational cost for analyzing HSM, we propose an efficient approximation scheme based on Importance Sampling and Empirical Interpolation Method. We illustrate our method using two examples - a Molecular Dynamics simulation for Krypton and a pharmacokinetic/pharmacodynamic model for cancer drug.


Introduction
Hierarchical Bayesian Model (HBM) is a powerful modeling tool extended from the classical Bayesian framework. It has been used to solve many difficult practical problems, such as modeling heterogeneous data, calibrating system with multiple objectives and inducing sparsity in a model. Its applications span across various fields, such as, social science, economics, physics, medical science, civil engineering, etc [28,9,2,14].
In a classical Bayesian setting, users define a stochastic model with parameters θ for a forward problem that predicts a quantity of interest and a prior distribution of θ. When data is available, the Bayes' Theorem is used to solve the inverse problem by finding the posterior distribution of θ. For a complex system in reality, we further parameterize the stochastic model and the prior with hyperparameters. This extra level of parameters provide extra flexibility to a model, but generally requires more data to reach a well-posed problem. HBM refers to model classes with such a hierarchy (multiple levels) of parameters. Here, we define two types of HBM commonly found in the literature: a Hierarchical Prior Model (HPM) that further parameterizes the prior, and a Hierarchical Stochastic Model (HSM) that further parameterizes the stochastic model (or known as the likelihood function when evaluated at a given data).
HPM is a well-studied subject in the machine learning and compressive sensing community in the context of Sparse Bayesian Learning (SBL). Its applications include imaging [23,6], gene selection [1], Bayesian compressive sensing [18,17], etc. It is an effective tool to solve ill-posed regression problems by adding extra constraints in the prior distribution under the Bayesian framework. The well-known Bayesian Optimization method can also be interpreted using this model structure. On the other hand, the origin of HSM can be traced back to the multilevel model in the statistics community. One classic application is to analyze test results collected from multiple schools [11]. [24] and [16] are early publications on Bayesian analysis for multilevel regression, and Congdon [8] summarizes recent developments on this topic. Since most of the computational demand in Bayesian analysis comes from evaluating the stochastic model, HSM has a significantly larger computational cost than HPM. Although recently, there is increasing amount of HSM applications, the diversity of HSM uses is not comparable to HPM. However, recent developments in parallel computing opens up new opportunities to apply HSM to more complicated problems. In this paper, we discuss important theoretical implications of HSM and demonstrate an efficient approximation scheme with practical applications.
Despite the frequent use of the terminology "HBM" in the literature, the distinction between what we called the HSM and the HPM is often omitted. The overlapped usage of "HBM" has even caused confusion. For example, Guha et al. [13] studied the HPM. Nevertheless, the authors used [12] as a reference, which the type of HBM mentioned in [12] is actually the HSM. Hence, in this paper, we begin with a comparison between HPM and HSM in Section 2. Then, we turn our focus back to HSM and study its theoretical implications based on polynomial regression in Section 3. To tackle the high computational cost for analyzing HSM in practice, we propose an efficient approximation based on the idea of Importance Sampling and Empirical Interpolation Method [3] in Section 4. Section 5 includes two realistic examples using our approximation method. Finally, we give some concluding remarks in Section 6.

Classification of HBM: HPM versus HSM
For the ease of illustration, we begin with defining the notations that will be repeatedly used in this paper. To predict a quantity of interest y, one may have a stochastic forward model F (x, θ, ) with model parameters θ and the stochasticity part . This function defines the distribution p(y|x, θ). One of the most common example of F (·) is: where y is chosen to follow a zero mean and σ y standard deviation Gaussian distribution, N (0, σ y ), for computational convenience. Given a set of input and output data pairs D = {(x i ,ŷ i )|i, . . . , N D }, the posterior distribution of θ is inferred by the Bayes' Theorem: where the likelihood p(D| θ) p(ŷ 1 , . . . ,ŷ N D |x 1 , . . . ,x N D , θ) is the probability of observing the data values D from the predictive model, the prior p( θ) is the initial belief of the values of θ, and the evidence (or marginal likelihood) p(D) p(ŷ 1 , . . . ,ŷ N D |x 1 , . . . ,x N D ) is a critical term used in model selection [5,4]. Then, a robust posterior prediction for an unobserved input-output pair (x 0 , y 0 ) can be obtained by: p(y 0 |x 0 , D) = p(y 0 |x 0 , θ)p( θ|D) d θ. ( HBM adds an extra level of hyperparameters ψ to this classical Bayesian model. This new hierarchy can be added to either the prior (HPM) or the stochastic model/likelihood (HSM), which will result in completely different model classes suitable for different problems. The major distinction of the two models comes from the different structure of information dependencies between all the stochastic variables. The graph representations shown in Figure 1 are very useful to display and study such relations [19]. The arrows denote the directions of the forward models, which implies that Bayes' Theorem is needed for inference in the opposite direction. Any nodes without an incoming arrow requires a prior distribution.
For HPM, an extra node is added to the starting nodes in order to parameterize the prior. This operation does not change the information dependency between the existing parameters. The posterior joint distribution of θ and ψ is: If we integrate out ψ in the analysis, this is equivalent to recovering the classical Bayesian setting from HPM by choosing p( θ) = p( θ| ψ)p( ψ) d ψ, because  HPM does not affect the stochastic model part. Hence, the hyperparameters in HPM behave as latent variables. In most applications of HPM, we perform optimization of ψ by maximizing p(D| ψ), instead of the robust treatment of ψ. For example, in the context of SBL, zero-mean Gaussian priors with different variances are chosen for the coefficients in a linear regression problem. This prior is called the Automatic Relevance Determination (ARD) prior that is proven to induce sparsity [20,26]. A commonly used algorithm in SBL, called the Relevance Vector Machine, is based on this optimization setup [25,27].
One can interpret this procedure as a Bayesian model selection problem in the continuous space. For HSM, the hyperparameters affect the basic structure of the stochastic model. As shown in Figure 1b, each prediction is defined by a separate set of model parameters, and these parameters are correlated through the hyperparameters ψ. Here, the original model parameters behave as latent variables, and ψ becomes the essential variables for future predictions y 0 . The posterior distribution of ψ is: Here, we denote the full data set as D = {D i |i = 1, . . . , N D }, where D i represents the observed data values for y i . The choice of p( θ i | ψ) is flexible. For example, a statistical model is chosen in [12] and [8] and ψ is called the hyperparameter vector. The statistical model can be seen as a common prior for all θ i . This explicit modeling for the prior separates different uncertainties by isolating out the uncertainty of θ i across multiple groups of predictions. It can give better predictions and also be related to causal inference [11]. Another example of such a hierarchical structure is found in a state-space model, where θ i is the state variable and p( θ i | ψ) represents the underlying stochastic process. Therefore, p( θ i | ψ) can also be a theoretically informed model, instead of a purely empirical or statistical model. In the current literatures, HPM is mainly used as a tool for optimal prior selection, while HSM is for explicit modeling of the uncertainty of model parameters across multiple groups of predictions for the data. However, HSM has the potential for more sophisticated modeling depending on the choice of p( θ i | ψ), though it comes with a significantly larger computational cost due to the extra integral shown in Equation 5. We note that since HPM and HSM are two exclusive types of HBM, in theory, they can co-exist in a complex HBM and may have more than one level of hyperparameters. It is the limit of computational power that constraints the usage of HBM to be only one level of hypereparameters. In this paper, we study the important theoretical implications of HSM and present an efficient and flexible approximation method for performing Bayesian analysis with HSM. We note that a similar comparison can be found in [22], but they present it purely from the perspective of a likelihood function of ψ, instead of the HBM.

Theoretical Implications of HSM
Different assumptions on the uncertainties in a system can be represented using different hierarchical and non-hierarchical models. In this section, we verify that the effect of Occam's Razor in the Bayesian model selection framework is applicable to both types of models. Moreover, we study the theoretical implications of using a HSM in three important aspects: 1. Separation of different types of uncertainties in a system (Section 3.1) 2. Identification of correlation between predictions or data (Section 3.1.2) 3. Uncertainty quantification of reduced order models (Section 3.2) For computational accuracy, we demonstrate our results using polynomial regression because of the availability of many analytical solutions during the Bayesian analyses.

Separating different types of uncertainties using HSM
In this study, we demonstrate how HSM can capture different types of uncertainties using models that are embedded with different sources of stochasticity. We generate synthetic data based on two uncertainty models: (1) zero-mean Gaussian error y added to the function f (x, θ), which can represent measurement noise in practice; (2) zero-mean Gaussian error θ added to the model parameters θ, which can represent inherent model uncertainty due to, for example, insufficient knowledge of theoretically informed models or environmental variations across multiple experiments. We use a linear function f (x, θ) = θx (θ and x are scalars) and all data points are generated with independent y and θ . Three types of synthetic data are considered based onθ, a fixed value of the model parameter θ: where N (z|µ, σ 2 ) denotes a Gaussian distribution for z with mean µ and variance σ 2 . We setθ = 1,σ θ = 0.5 andσ y = 0.2. Each data set contains 1000 data points independently generated from a uniformly distributed x value between preset bounds and random errors y and θ from the corresponding Gaussian distributions. We generate all three types of data twice, once with x between 0 and 1 and once with x between 0.4 and 1 (see Figures 2 and 3, respectively).  Given the three types of data, we perform Bayesian model class selection on four different stochastic model classes: 1. Non-hierarchical model, M 1a : This model class includes a uniform prior for θ between -1 and 3, and a Gaussian likelihood for θ, i.e., p({x i , y i }|θ, σ y ) = N (y i |f (x i , θ), σ 2 y ) for any i. Here, σ y is also treated as a parameter to be inferred. We use a uniform prior for σ y between 0.001 and 1.
3. Zero noise HSM, M 2a : This is an HSM with p(θ i | ψ) modeled as the same Gaussian prior used for θ in M 1b for all i. No additive error is assumed in this model. This implies that the likelihood function of θ i is a delta function centered at the given data, i.e., The priors of the hyperparameters are the same as in M 1b .
4. Full HSM, M 2b : This is the same HSM as M 2a except that a Gaussian additive error is assumed in this model. Hence, the likelihood function . We use the same uniform priors as in M 1b for σ y , µ θ and σ θ .
Appendix A includes all derivations of the analytical expressions used for estimating the evidence p(D|M k ), the marginal posterior probability density function (PDF) for the model parameters p(θ|D, M k ), and the posterior PDF for the hyperparameters p(ψ|D, M k ), for a given model class M k . We use Monte Carlo Simulation to draw posterior samples of the hyperparameters for the "naive" estimation of the evidence. 10K samples are used to obtain accurate estimates with small variances. Table 1 and 2 summarize the results of the Bayesian inference and model selection for all six data sets (three with x between 0 and 1, and three with x between 0.4 and 1). The Bayesian model selection framework is capable of selecting the corresponding model used to generate the data. All models have a relatively good estimate forθ given by E[θ|D]. The major difference appears in the accuracy of estimatingσ θ andσ y given by Std[θ|D] and E[σ y |D], respectively. In general, in terms of accurately estimatingσ θ andσ y , M 1a and M 1b always prefer putting uncertainty into σ y and thus perform well for D 1 only. M 2a cannot handle additive error and thus performs well for D 2a only. M 2b is the most flexible model and it performs relatively well for all data sets D 1 , D 2a and D 2b .

Results and discussion
We note that the results of M 1a and M 1b are extremely similar. This is because their only difference is the prior of θ, as explained in Section 2. In this case, the marginalized prior p(θ) in M 1b is slightly larger than the prior in M 1a . Hence, the evidence of both models are almost the same in all cases.
Furthermore, we observe that in the case of D 2b with x between 0.4 and 1, M 2a and M 2b have a similar posterior model probability. Intuitively, the data with small x values are more sensitive to additive noise because the noise-tosignal ratio is much higher than in the case of larger x values. When such data is not available, it is challenging to distinguish between D 2a and D 2b . This can be proved visually by observing the similarity between D 2a and D 2b in Figure 3 as compared to the one in Figure 2.
Our results suggest using M 2b for model calibration, when no prior knowledge indicates that additive error is irrelevant. This HSM can appropriately separates the two types of uncertainties in our case study into the additive and non-additive parts based on the observed data. If computational power allows, Bayesian model class selection should always be performed among a set of different candidate models, as illustrated in this study.

Effect of grouping
In the previous section, all synthetic data is generated independently based on different stochastic models. In practice, we may group some of the predictions to be modeled under one set of parameter values θ i , where i is the index for each group of predictions. For example, if an experiment is repeated in different laboratories, we may model the data from the same laboratory using the same set of model parameters, assuming the data shares the same environmental factors during the experiment. When the prediction grouping is not known, we may want to know the most plausible grouping, which defines a unique HSM.
We use the same linear function f (x, θ) = θx as before, but generate a new set of data D. First of all, five random samples θ (i) for i = 1, . . . , 5 are drawn from N (θ|μ θ ,σ 2 θ ). Then, 11 data points are generated for each θ (i) to form a data set D i by using the stochastic forward model y = f (x, θ (i) ) + y , where y is random error drawn from N ( y |0,σ 2 y ) independently for each data point. In 1.000 0.500 Figure 4a shows the five data sets used in the section. The sample mean and standard deviation of θ (i) are 1.0191 and 0.4079, respectively.
To study the effect of grouping predictions in the HSM, we set up six different HSMs based on different information dependency between the predictions. We perform Bayesian model class selection to find the most plausible HSM to describe the data. The evidence is calculated similarly to M 2b in Section 3.1. Details of the calculation can be found in Appendix A. The six data groupings are: 1. Actual grouping, M 1 : group predictions according to Figure 4a.
2. Constant x grouping, M 2 : predictions with the same x value are grouped together. This case represents grouping data according to the wrong variable. Figure 4b shows the resulting grouping.
3. 1/2 error grouping, M 3 : starting from M 1 , the predictions for each data set are further divided into 2 sets by the center line of the data points in the set. This center line is defined by the mean of θ inferred from each data point in the set. Figure 4c shows the resulting data grouping. 5. Single prediction, M 5 : each prediction is treated as an independent data set (same as Section 3.1.1).
6. Random grouping, M 6 : predictions for the data points are randomly collected into five groups. Figure 4d shows the resulting grouping.   Table 3 summarizes the results of the Bayesian inference and model class selection for all six HSM. All models have a good estimate forμ θ (corresponds to E[θ|D]). However, M 2 and M 6 , representing two completely wrong groupings, have significantly lower estimates forσ θ (given by Std[θ|D]), but larger estimates for the mean ofσ y (given by E[σ y |D]). Hence, they have a very low log-evidence value. M 1 has the closest estimates ofμ θ ,σ θ andσ y , but it is not the most probable model class among the six. M 3 is the most probable model class even though it has a low value of the estimates for bothσ θ andσ y . We note that when the total number of data points is fixed, the more groups there are, the smaller the average number of predictions in a group is. This affects the likelihood of ψ, which equals the product of p(D i | ψ) for each data set D i (see Equation 5). When the number of p(D i | ψ) increases, the value of each evidence term may decrease. This decrease is due to a less peaked p(D i | θ i ) as the number of data points in the set decreases. The final evidence p(D|M k ) for a given data grouping M k is a tradeoff between these two factors. Starting from M 1 being modified to M 3 , if we repeat the process many times, eventually we will reach M 5 . Therefore, M 1 , M 3 , M 4 and M 5 can be considered as a sequence of similar data grouping methods. In the end, the tradeoff suggests that M 3 is the most probable grouping. This implies that model class selection based on the evidence may not necessarily lead to the actual grouping, if it exists. Overall, it tries to minimize uncertainties in all parameters. On the other hand, if we choose the actual grouping, the uncertainty quantification will be accurate.

Uncertainty quantification for reduced order models using HSM
A reduced order model simplifies the information obtained in the original model. This loss of information can be treated as a source of uncertainty in the Bayesian framework. The structure of this type of uncertainty can be defined if the mapping between the original model and the reduced order model is completely known. In most cases, however, the mapping is not well-defined or too complicated to be studied directly. Therefore, there is nono clear way to model such kind of loss of information. In this section, we compare the performance of using different model classes in Section 3.1 to represent the reduced order model  uncertainty. This study is based on fitting second and third degree polynomial data with linear functions.

Problem setup for polynomials
We consider a total of 16 sets of data: 8 sets from a quadratic function f (x) = x 2 and 8 sets from a cubic function f (x) = x 3 . For each type of function, we generate 20, 50, 100 and 200 data points twice (with and without noise), i.e., a total of 8 data sets. The x values of the data points are generated randomly in the interval [−1, 1]. Additive Gaussian error with standard deviationσ y = 0.1 is chosen for the noise. Figure 5 shows all 16 data sets used in this study. We perform Bayesian model class selection and posterior robust prediction using three models described in Section 3.1: M 1a (denoted as M 1 in this section), M 2a and M 2b . Appendix A includes all derivations for the analytical expressions used to estimate the log-evidence values and posterior robust predictions p(y|D, M k ) for a given model M k . We evaluate the posterior probability of observing each grid point (x,ŷ) on a 2D fine grid of (x, y) using the analytical expressions in order to construct the distribution of the posterior robust prediction. Table 4 summarizes the results of model class selection for this study. M 2b is the most probable model class for data with additive noise and M 2a is the most probable model class for data without any additive noise. M 1 is the least significant model in most cases because the error caused by the reduced order model is very non-linear. We do not expect that the additive noise alone is sufficient to explain the error. In this study, the sample size (number of data points) does not play a significant role in the results. Figure 6 shows one set of the prediction results (the results are not sensitive  These results suggest that the HSM with additive error is preferred as the first test for uncertainty quantification of reduced order model. The flexibility of such a model class to handle two types of uncertainty simultaneously (additive error and embedded error in the model parameters) results in a higher chance of discovering the underlying uncertainty structure of a reduced order model.

Efficient Approximation of HSM
The posterior distribution p( ψ|D) plays an important role in the HSM. Equation 5 shows that the calculation of the likelihood p(D| ψ) involves evaluations of multiple integrals, which correspond to the evidences for each data set D i conditional on the hyperparameters ψ. This leads to an extremely large computational cost for sampling from p( ψ|D), as well as estimating the evidence of the model, p(D). Current approaches include using conjugate pairs for analytical results [8], approximating the integrals with Laplace Asymptotic Approximation [29], or using some advanced Markov Chain Monte Carlo techniques [21]. These methods are still restrictive for either studying many complex systems or the efficiency is not very scalable to handle extra data sets.
In this section, we present an efficient approximation method based on a special use of Importance Sampling. Most current research efforts focus on  Figure 6: Posterior robust prediction for different models in different cases (sample size = 100). The purple dash lines denote the mean prediction and the grey area encloses 90% of the total probability density for the predicted value. Blue crosses are the data points. The four columns correspond to those in Table 4. the HSM that adds only one extra level to the classical Bayesian model. Our method can be intuitively and efficiently extended to any complex HSMs with more levels of parameters. Moreover, we can lay out a standard procedure for fully analyzing any HSMs based on this method.
For ease of illustration, we demonstrate our method based on the commonly used stochastic model shown in Equation 1. In practice, there may be more complicated HSM than the one shown in Figure 1b. Here, we consider two possible alternatives (see Figure 7): (1) independent additive error parametersapplicable to heterogeneous data with different likelihoods [30], and (2) common additive error parameters -applicable to data sets that are expected to have the same measurement errors [29].
We introduce our method by first applying it to the basic HSM shown in Figure 1b, which represents the cases that either σ y is known or it is uncertain and included with the other uncertain parameters θ. Then, we extend the idea to the two alternative HSMs shown in Figure 7.

Basic HSM
This model class considers σ y to have a known value or to be treated as part of θ. Hence, we do not infer the distribution of σ y from the data D separately. Figure 1b shows the graphical representation for this case. In order to separate it from the two alternative HSMs in Figure 7, we denote this model class as M HS1 . The posterior distribution, the likelihood and the evidence for ψ are: Our idea is to approximate the integrals by Importance Sampling with a special choice of the proposal PDF that can significantly reduce the total computational cost. In many cases, the likelihood p(D i |θ i , M HS1 ), which is related to the evaluation of f (x, θ), dominates the computational effort. Our method begins by performing Bayesian inference for each data set D i using the same likelihood p(D i |θ i , M HS1 ), but a different prior (choice of such prior is discussed later). To be more specific, we draw samples { θ Then, we can approximate p(D i |ψ, M HS1 ) based on IS with proposal distribu- As a result, we only need to perform classical Bayesian inference once for each data set D i (draw posterior samples and estimate the evidence p(D i |M i )). Then, the hierarchical analysis comes as a post-processing of the results. A very good tool for such a Bayesian inference is the Transitional Markov Chain Monte Carlo (TMCMC) method, where the evidence comes as a by-product of drawing posterior samples and the algorithm is inherently parallel [7]. The likelihood p(D|ψ, M HS1 ) can then be estimated by: In other words, the evidence ratio between M HS1 with given hyperparameters and M i can be approximated by the mean of the prior ratio between the two models over the posterior θ i samples of M i . This approximation suffers the same problem as IS, i.e., the variance of the estimate depends strongly on the closeness of the integrand and the proposal distribution. In our case, they are exactly the same if p( θ i |M i ) = p( θ i | ψ, M HS1 ). Therefore, we should choose p( θ i |M i ) to minimize its difference to p( θ i | ψ, M HS1 ) for all ψ weighted by the prior p( ψ|M HS1 ). If we use the Kullback-Leibler divergence D KL (p||q) as a measure of difference between two distribution p(x) and q(x), this implies the choice: A large number of posterior samples from p( θ i |D i , M i ) can also reduce the variance of the IS estimate. If this is computationally feasible, an alternative for p(θ i |D i , M i ) would be a uniform distribution that covers the significant regions of p( θ i |M HS1 ) = p( θ i | ψ, M HS1 )p( ψ|M HS1 ) d ψ. We note that the evidence will be insignificantly small for ψ values that lead to the prior p( θ i | ψ, M HS1 ) having a mode well away from the maximum of the likelihood p(D i | θ i , M HS1 ). Hence, accuracy is not important in this case. In the opposite case, IS is a good approximation. The inaccuracy problem is important only for the middle case, i.e., the prior mode is not too far and not too close to the maximum of the likelihood. The range of ψ leading to this important case depends on the posterior sample size N s,i for each data set D i .
We note that this approximation method provides information about Bayesian inference for each individual data set first. This information is useful to compare with the HSM analysis to give further insight about the system of interest. Furthermore, our method is very efficient for introducing extra data sets because we do not need to rerun the whole problem. Instead, the overhead is only the classical Bayesian inferences for the extra data sets and a post-processing step that is not computationally intensive.

HSM with independent additive error parameters
This model class considers σ y for the predictions of each data set to be independent of each other. Figure 7a shows the graphical representation for this case.
p(D| ψ, σ y,1 , . . . , σ y, The idea for M HS1 does not apply to this model class directly because both the likelihood p(D i | θ i , σ y,i , M HS2 ) and the prior p( θ i | ψ, M HS2 ) will be affected by the sampling of ψ, σ y,1 , . . . , σ y,N D . However, if prediction is the ultimate goal, Figure 7a implies that σ y,1 , . . . , σ y,N D are not relevant as long as the posterior distribution of ψ is known. Then, we can again use IS to estimate p( ψ|D, M HS2 ) directly.
We note that this is a special case of M HS1 with p( θ i , σ y,i | ψ) = p( θ i | ψ)p( σ y,i ), i.e., σ y,i is independent of ψ. Similar to M HS1 , we first draw samples {( θ Following the procedure for M HS1 , we can derive that: Again, the HSM analysis comes as a post-processing step and p(D| ψ, M HS2 ) = As a result, we can directly obtain posterior samples from p( ψ|D, M HS2 ) using MCS or MCMC methods combined with Equation 19. We note that the estimation can be further simplified and more accurate when choosing p(σ y,i |M i ) = p(σ y,i |M HS2 ) for all i.

HSM with common additive error parameters
This model class considers a single value of σ y to be shared by all predictions. Figure 7b shows the graphical representation for this case. We denote this model class as M HS3 . In this model class, σ y affects future prediction as well. Hence, the approach for M HS2 will not work. We adopt the idea from EIM to solve the problem. The basic concept of our approach is to approximate the likelihood function p(D i | θ i , σ y , M HS3 ) as a linear sum of multiple copies of the same term, but with fixed σ y values. Then, we can perform the classical Bayesian analyses for those fixed σ y cases and use the same IS approach as in M HS1 and M HS2 to approximate the HSM analysis for M HS3 . First, we look at how to use the EIM idea to estimate the posterior distributions in M HS3 . Then, we discuss how to train the linear approximations for the likelihood p(D i | θ i , σ y , M HS3 ).
We denote g l ( θ i ) = p(D i | θ i , σ y,l , M HS3 ) as a basis function for some fixed value of σ y,l for l = 1, . . . , L. EIM assumes the following approximation for the fixed number of bases L and some coefficients α l ( σ y ): Then, we can approximate the likelihood p(D i | ψ, σ y , M HS3 ) = p(D i | θ i , σ y , M HS3 )p( θ i | ψ, M HS3 ) d θ i using Equation 20: Since σ y,l is fixed, we can apply the same IS approach as in M HS2 and M HS1 to the integrals in Equation 21. We draw samples { θ (j) Then, we use the proposal q l,i ( θ i ) = p( θ i |D i , σ y,l , M i ): As a result, the likelihood p(D| ψ, σ y , M HS3 ) can, then, be estimated by: With this analytical expression for the approximation, the remaining problem is how to pick σ y,l and how to find α l ( σ y ) when estimating the posterior distribution p( ψ, σ y |D, M HS3 ).

Training basis functions
Hesthaven et al. [15] suggest using an adaptive greedy algorithm to select the set of basis functions g l (θ i ). Based on this idea, we develop an algorithm that simultaneously selects basis functions and collects posterior samples used in Equation 24. The greedy algorithm reduces the maximum absolute value of the error term e l ( θ i , σ y,l ) to some specified threshold˜ lim over training sets of θ i and σ y , denote as Θ i and Σ y , respectively. The efficiency and accuracy of the algorithm on error estimation is a tradeoff that depends on the size of the training sets. The error term corresponding to L basis functions is defined as: Starting with some initial training sets Θ initial i and Σ initial y and an initial set of basis functions G initial , a new basis function with a corresponding σ y value is chosen from Σ y to maximize the error term e L . Also, we record the θ i values that maximize e L for the chosen σ y . These values will be used in the "online" estimation stage discussed later. Then, a set of posterior samples { θ i } L is drawn using any MCMC methods with the chosen σ y , and the samples are added to the current training sets Θ i . This process is repeated until the maximum error given Θ i is below the threshold e lim . In the end, we obtain a set of bases defined by ( θ i,l , σ y,l ), l = 1, . . . , L, the posterior samples { θ i } l and the evidence values p(D i | σ y,l , M HS3 ) for estimation of the approximate hierarchical Bayesian inference. Algorithm 1 summarizes this new adaptive training procedure.

Algorithm 1 : Adaptive greedy algorithm
Obtain initial σ y,1 , . . . , σ y,L0 for basis functions g 1 , . . . , g L0 from G initial Obtain initial θ i,1 , . . . , θ i,L0 for the corresponding basis functions from G initial Initialize counter l ← L 0 Initialize training sets Θ i = Θ initial i and Σ y = Σ initial ỹ L0 = max σy∈Σy max θi∈Θi e L0 ( θ i , σ y ) while˜ l >˜ lim do l ← l + 1 σ y,l = argmax σy∈Σy max θi∈Θi e l−1 ( θ i , σ y ) By enriching the training set Θ i , the error estimate e L ( θ i , σ y ) becomes more accurate and thus improves the EIM approximation. In many cases, most of the computational time for calculating p(D i | θ i , σ y , M HS3 ) is spent on evaluating the function f (x, θ i ). Because this evaluation is already done during MCMC sampling, the overhead of using a larger training set Θ i is relatively small. However, it is important to note that the improvement of the error estimate may saturate. This occurs when the different sets of posterior samples for different σ y values all have the same high probability regions. In this case, there will be many redundant samples in the training set clustering around the peaks of p(D i | θ i , σ y , M HS3 ). To further improve the efficiency of this training algorithm, we can control the expansion of Θ i such that we only add samples that are significantly different from the samples in the current Θ i . A simple implementation is to monitor the spread of the chosen σ y because similar σ y values represents a similar likelihood function value p(D i | θ i , σ y , M HS3 ) and thus similar posterior samples are expected. Using the same argument, we suggest constructing the initial sets Θ initial i and G initial based on extreme values of σ y (e.g., maximum and minimum values of σ y in 1D case). Algorithm 2 constructs the initial sets based on this suggestion.
Pick a fine grid for Σ initial y Initialize σ y,1 , . . . , σ y,K as a sequence of extreme σ y values ordered in ascending "peakness" of p(D i | θ i , σ y , M HS3 ), i.e., more probability content concentrated in a small region of θ i values (only qualitatively, need not to be accurate). This is the initial set of σ y values that define the initial basis functions g 1 ( θ i ), . . . , g K ( θ i ) in G initial Initialize Θ Initial i ← empty set or a set of sparse grid points

Online estimation
Once we obtain the set of bases ( θ i,l , σ y,l ) that defines the basis functions g l ( θ i ), where l = 1, . . . , L and L is the total number of bases, we are ready to estimate p(D i | ψ (s) , σ Once α l is solved, we can include the recorded posterior samples { θ i } l and evidence values p(D i | σ y,l , M i ) corresponding to basis g l ( θ i ) to perform fast estimation of p(D i | ψ, σ y , M HS3 ).

Numerical issue
In Equation 26, note that the vector {P n } and each column of the matrix [g nl ] represents the likelihood value p(D i | θ i,n , σ y , M HS3 ) for different values of σ y . as: If we choose c α l = c P /c g l for all l = 1, . . . , L, the inversion problem becomes: We pick c P and c g l to be the maximum value of the corresponding column vector. Then, the numerical issue due to scaling may not occur in Equation 28. We can first solve the inversion problem in Equation 28 and recover the coefficients α l by {α l } = {c α 1α1 · · · c α Lα L } T where c α l = c P /c g l . Also, it may be useful to store the values in log-scale during the entire computation process.

Illustrative Examples
We test our method using two examples: molecular dynamics simulation of Krypton and pharmacokinetic/pharmacodynamic (PKPD) modeling of a cancer drug. Both examples are taken from previous studies. Their results are used to validate the efficiency and accuracy of our method.

Molecular Dynamics: Krypton
Wu et al. [29] proposed using the HSM (called the HBM in the paper) to calibrate the Lennard-Jones (LJ) potential ( LJ and σ LJ ) of Krypton based on viscosity data D reported from nine different laboratories. In that paper, the authors approximated the integrals in Equation 7 by Laplace Asymptotic Approximation (LAA). Since the joint posterior distribution of LJ and σ LJ for each data set is close to a Gaussian distribution, LAA is a good approximation for this particular case. We re-visit the problem with the same data and the same model class setup as described in [29]. Instead of the LAA approximation, we use our proposed method to estimate the joint posterior distribution of LJ and σ LJ and the posterior robust prediction.
Following the assumptions in [29], we use the HSM with common additive error. First, we build the EIM basis functions based on Section 4.3.1. We begin with a training set of LJ and σ LJ from a coarse grid of 256 points within the boundaries 100 ≤ LJ ≤ 400 and 0.2 ≤ σ LJ ≤ 0.5. A fine grid of 172 points for σ y is used between 0.0001 and 0.5, uniformly distributed in log scale. We obtain a total of 50 basis functions for each data set to achieve an error less than 0.001%. Figure 8 shows the results of the EIM approximation. We observe that the bases of LJ and σ LJ coincide with the high likelihood values in the domain. This is consistent with our intuition that important regions of the likelihood should be included in the online estimation for better accuracy.  For each EIM basis, 2500 posterior samples of LJ and σ LJ are recorded to be used in the post-processing step of our HSM analysis. We use BASIS to draw the samples because it is highly parallel and the evidence is a by-product of the algorithm [31]. Then, we draw the posterior samples of the hyperparameters ψ with the post-processing step in our method. BASIS is once again used to draw 1000 posterior samples of ψ. The posterior distribution of LJ and σ LJ is approximated by N posterior samples of ψ: Figure 9 shows the results of the HSM analysis. We observe that our results are consistent with the results reported in [29]. This provides a verification of the accuracy of our method.  calibrated using two types of clinical data: the plasma concentrations of free aflibercept and bound aflibercept. The data can be separated into 6 groups, each corresponding to a different dosage of the drug. We apply the HSM structure in Figure 1b to the basic model used in [10], which has three PK model parameters and five prediction error parameters. The prediction of the two types of data using these parameters are calculated by first running the system to reach steady state, and then imposing the corresponding drug dosages. The three model parameters {k EC , k N , k T } represent the secretion rate of VEGF in the blood, normal tissue and tumor compartment, respectively. The five prediction error parameters {σ SS , σ BAs , σ BAf , σ F As , σ F Af } represent standard deviations of zero-mean Gaussian distributions for five different types of prediction errors: σ SS -error on the steady state prediction, σ BAs -error on the bound aflibercept prediction scaled by time, σ BAf -error on the bound aflibercept prediction without any scaling, σ F As -error on the free aflibercept prediction scaled by time, σ F Af -error on the free aflibercept prediction without any scaling. We choose independent Gaussian prior distributions for all of the eight parameters in log-10 scale. Figure 10 and Table 5 compare the posterior distribution of the classical Bayesian model used in [10] with the posterior distribution of the basic HSM presented in this paper. We observe that fitting all data at once (the classical Bayesian model) leads to a significantly larger posterior variance for the model parameters, while the HSM leads to a larger posterior variance for the prediction error parameters. This implies that the knowledge about the model parameters is more transferable across different dosages than the knowledge about the prediction error parameters. We note that the steady state prediction is not affected by the change of dosage. Indeed, we observe that the inferred value of σ SS in the HSM is similar to the value for the classical Bayesian model. The slight increase of the CV of σ SS for the HSM case is expected because the data used for inference in the HSM is divided into groups of smaller data sets. Moreover, the scaled standard deviations σ BAs and σ F As are larger, and the fixed standard deviations σ BAf and σ F Af are smaller in the HSM. This is also ex-pected because in the classical Bayesian model, predictions for different dosages may have different level of sensitivity to the errors scaled by time. Hence, the classical Bayesian model tends to have a larger σ BAf and σ F Af than σ BAs and σ F As . As a result, the posterior robust prediction based on the HSM ( Figure  11) has completely different prediction errors for both bound aflibercept and free aflibercept compared to the one in [10].

Conclusion
We demonstrate the benefits of using the hierarchical Bayesian framework in complex systems. Because of the double usage of the term HBM in the literature, we first explain the distinction between two very different HBMs: the HPM and the HSM. Then, we focus on studying the HSM, which has many interesting theoretical implications, as well as many computational challenges in practice.
Based on examples of polynomial functions, we suggest that the HSM is capable of explicitly separating different types of uncertainties in a system, and can be an effective tool for modeling uncertainty of reduced order models. In order to apply HSM to practical problems, we propose an efficient approximation method to tackle the high computationally cost associated with using HSMs. Our method is a "bottom-up" approach that begins by drawing posterior samples of the model parameters for each data set. Then, we use a post-processing step to   move up the hierarchy of the parameters for drawing posterior samples. Building on the basic HSM, we demonstrate our method using two alternative HSMs that have a more complicated hierarchical structure. Lastly, we validate our method using two illustrative examples based on previous studies: a molecular dynamics simulation of Krypton and a PKPD model of a cancer drug. The HSM is a convenient and effective tool to build the stochastic model/likelihood for a complicated system. It also opens up new types of analyses and it results in different conclusions as compared to classical Bayesian inference. Our approximation method provides a standard procedure for analyzing hierarchical models. Beginning with an analysis for each individual data set, our method allows us to move up the hierarchy of the parameters efficiently to potentially extract more in-depth information about the system of interest. Moreover, it is very efficient for sequentially received data sets because our estimation scheme is a computationally fast post-processing step. So far, we have mainly used a purely statistical (empirical) model for p( θ| ψ). In future work, we plan to employ other choices for modeling p( θ| ψ), which may motivate the development of new modeling tools.

A Derivations
We consider a set of data D structured as D = {D i |i = 1, . . . , N D } with each subset of data D i = {(x j,i , y j,i )|j = 1, . . . , N Di }, where (x j,i , y j,i ) denotes a single data point. Using vector notations, we denote Hence, x and y are vectors with N d elements.
Starting with a general formulation, we denote θ as a vector of all model parameters and ψ as a vector of all hyperparameters. Repeatedly using Bayes' Theorem and the total probability theorem, we can derive the following expressions for a given model M k : We note that usually p(M k ) is constant for all k because we do not want to introduce bias to any model before data is available. Hence, p(M k |D) ∝ p(D|M k ). Also, θ is simply a scalar θ in our study and the likelihood standard deviation σ y is inferred separately. In the following, we derive the analytical expressions for the posterior distributions, the evidences and the robust-posterior predictions of the different models.

A.1 Non-hierarchical model, M 1a
In this model, all data points are independent when θ and σ y are known. By assuming a Gaussian distribution for the likelihood and a uniform prior U (θ) = 1/c θ , we can derive that: Using Bayes' theorem, the posterior distribution is proportional to only the likelihood when the prior is uniform. It can be obtained from completing square for the exponential part of the Gaussian likelihood: (38) Hence, the posterior distribution is also Gaussian: and the evidence term p(D|σ y , M 1a ) calculated using Equation 35 (substitute ψ by σ y ) is: Then, we use N s samples from Monte Carlo Simulation to estimate the posterior of σ y and the model evidence: The marginalized posterior of θ is estimated by substituting Equation 39 and 41 into Equation 31, which we substitute ψ by σ y : As a result, the statistics of the marginalized posterior of θ and σ y can also be estimated (E[·|D] -posterior mean; Std[·|D] -posterior standard deviation) based on the MCS samples: To perform robust prediction of a new point (x,ŷ), we need to evaluate p((x,ŷ)|D, M 1a ). Again, we can use the posterior samples from Equation 41: where p((x,ŷ)|σ As explained in Section ??, this model is essentially the same as M 1a , but with a different prior of θ. By assuming a Gaussian distribution for the likelihood and the prior of θ given the hyperparameters ψ = {µ θ , σ θ }, we can derive that: Using Bayes' theorem, the posterior distribution is proportional to the product of the likelihood and the prior. It can be obtained from completing square for the exponential part of the Gaussian distributions: Hence, the posterior distribution is also Gaussian: and the evidence term p(D|µ θ , σ θ , σ y , M 1b ) calculated using Equation 35 is: We use N s samples from Monte Carlo Simulation to estimate the posterior of µ θ , σ θ and σ y : where (µ Ns j=1 w j = 1 and the model evidence: The marginalized posterior of θ can be estimated by substituting Equation 50 and 52 into Equation 31: As a result, the statistics of the marginalized posterior of θ and σ y can also be estimated based on the MCS samples: In this model, we assume the same Gaussian prior of θ i for all i = 1, . . . , N D , where θ i corresponds to the data set D i ∈ D. The hyperparameters ψ = {µ θ , σ θ } define the mean and the standard deviation of the Gaussian prior, respectively. By assuming a delta function for the likelihood, we can derive that for each data set D i : i.e., p(θ i |D i , µ θ , σ θ , M 2a ) = 1, if there exists θ i such that y j,i = θ i x j,i for all j = 1, . . . , N Di 0, otherwise (58) And when the posterior equals 1, the evidence term p(D i |µ θ , σ θ , M 2a ) calculated using Equation 35 is: Otherwise, the evidence is zero. Note that the extra product of inverse of x j,i comes from the Jacobian of the delta likelihood function. Similar to M 1b , we use N s samples from Monte Carlo Simulation to estimate the posterior of µ θ and σ θ : where (µ (60) and the model evidence: In the HSFM, θ for future predictions (denoted as θ 0 in Section ??), is independent of the data when ψ is given. Hence, the marginalized posterior of θ depends only on the posterior of µ θ and σ θ and the prior of θ given µ θ and σ θ : As a result, the statistics of the marginalized posterior of θ can also be estimated based on the MCS samples: Note that an unique feature of the HSFM is that future predictions do not depend on the data once all hyperparameters are given. Hence, p((x,ŷ)|µ

A.4 Full HSFM, M 2b
In this model, we use the same setup as in M 2a for the hyperparameters ψ. By assuming a Gaussian distribution for the likelihood and the priors of θ i for all i = 1, . . . , N D conditional on the hyperparameters, we can obtain analytical expressions similar to the one in M 1b . One difference is that in this model, we now work with multiple data sets D i ∈ D, i = 1, . . . , N D . Following a similar derivation as in M 1b , the posterior distribution p(θ i |D i , µ θ , σ θ , σ y , M 2b ) is Gaussian: p(θ i |D i , µ θ , σ θ , σ y , M 2b ) = N (θ i |μ θ,i , (σ θ,i ) 2 ) and the evidence term p(D i |µ θ , σ θ , σ y , M 2b ) calculated using Equation 35 is: Similar to M 1b , we use N s samples from Monte Carlo Simulation to estimate the posterior of µ θ , σ θ and σ y : p(µ θ , σ θ , σ y |D, M 2b ) ∝ p(D|µ θ , σ θ , σ y , M 2b )p(µ θ , σ θ , σ y |M 2b ) p(D i |µ θ , σ θ , σ y , M 2b )p(µ θ , σ θ , σ y |M 2b ) where (µ (j) Ns j=1 w j = 1 and the model evidence: Similar to M 2a , the marginalized posterior of θ can be estimated by: As a result, the statistics of the marginalized posterior of θ and σ y can also be estimated based on the MCS samples: