Second-Order MaxEnt Predictive Modelling Methodology. I: Deterministically Incorporated Computational Model (2nd-BERRU-PMD)

Dan Gabriel Cacuci

Center for Nuclear Science and Energy, Department of Mechanical Engineering, University of South Carolina, Columbia, SC, USA.

DOI: 10.4236/ajcm.2023.132013                

Abstract

This work presents a comprehensive second-order predictive modeling (PM) methodology designated by the acronym 2^nd-BERRU-PMD. The attribute “2^nd” indicates that this methodology incorporates second-order uncertainties (means and covariances) and second-order sensitivities of computed model responses to model parameters. The acronym BERRU stands for “Best- Estimate Results with Reduced Uncertainties” and the last letter (“D”) in the acronym indicates “deterministic,” referring to the deterministic inclusion of the computational model responses. The 2^nd-BERRU-PMD methodology is fundamentally based on the maximum entropy (MaxEnt) principle. This principle is in contradistinction to the fundamental principle that underlies the extant data assimilation and/or adjustment procedures which minimize in a least-square sense a subjective user-defined functional which is meant to represent the discrepancies between measured and computed model responses. It is shown that the 2^nd-BERRU-PMD methodology generalizes and extends current data assimilation and/or data adjustment procedures while overcoming the fundamental limitations of these procedures. In the accompanying work (Part II), the alternative framework for developing the “second- order MaxEnt predictive modelling methodology” is presented by incorporating probabilistically (as opposed to “deterministically”) the computed model responses.

Keywords

Second-Order Predictive Modeling, Data Assimilation, Data Adjustment, Uncertainty Quantification, Reduced Predicted Uncertainties

1. Introduction

The results of measurements and computations are never perfectly accurate. On the one hand, results of measurements inevitably reflect the influence of experimental errors, imperfect instruments, or imperfectly known calibration standards. Around any reported experimental value, therefore, there always exists a range of values that may also be plausibly representative of the true but unknown value of the measured quantity. On the other hand, computations are afflicted by errors stemming from numerical procedures, uncertain model parameters, boundary and initial conditions, and/or imperfectly known physical processes or problem geometry. Therefore, knowing just the nominal values experimentally measured or computed quantities is insufficient for applications. The quantitative uncertainties accompanying measurements and computations are also needed, along with the respective nominal values. The discrepancies between experimental and computational results provide the basic motivation for performing quantitative model verification (meaning: “is the mathematical model solved correctly?”) and model validation (meaning: “does the model represent reality?”), which are essential components of “predictive modeling.”

Predictive modeling commences by identifying and characterizing the uncertainties involved in every step in the sequence of the numerical simulation processes that ultimately lead to a prediction. Predictive modeling comprises three key elements, namely model calibration, model extrapolation, and estimation of the validation domain. Model calibration addresses the combination of experimental and computational data and their uncertainties for the purpose of obtaining “best estimate” values for model parameters (to be used for updating the model’s parameters) and predicted model results, along with “best estimate” uncertainties (covariance/correlation matrices) for these “best-estimate” predicted parameters and results requires. Such a combination of computational and experimental information requires reasoning from incomplete, error-afflicted, and occasionally discrepant information, which includes: 1) errors and uncertainties in the data used in the simulation (e.g., input data, model parameters, initial conditions, boundary conditions, sources and forcing functions); 2) numerical discretization errors; and 3) uncertainties in (e.g., lack of knowledge of) the processes being modeled.

Under ideal circumstances, the result of predictive modeling is a probabilistic description of possible future outcomes based on all recognized errors and uncertainties. This probabilistic description enables the subsequent activity of “model extrapolation,” which aims at quantifying the uncertainties in predictions under new environments or conditions, including both untested regions of the parameter space and higher levels of system complexity in the validation hierarchy. The quantification of the validation domain underlying the models of interest requires estimation of contours of constant uncertainty in the high-dimensional space that characterizes the application of interest, including the identification of areas where the predictive estimation of uncertainty meets specified requirements for the performance, reliability, or safety of the system of interest.

The earliest activities aimed at extracting best-estimate values for model parameters by combining computational and experimental information using variational methods were initiated in the 1960s in the atmospheric and oceanographic sciences [1] [2] [3] and, in parallel, in the nuclear energy field [4] . In the earth, atmospheric and oceanographic sciences these activities were carried out under the name of “data assimilation” (DA) as described in well-known books [5] [6] [7] [8] , while in the nuclear sciences these activities reached conceptual maturity under the name of “nuclear data or cross section adjustment” [9] [10] [11] [12] . The “data adjustment” and “data assimilation” methodologies can thus be considered to have been the earliest systematic methodologies that embody the principles of “predictive modeling.

The fundamental criterion used in data adjustment and data assimilation methods is the least squares criterion, which is employed in a variety of deterministic (variational) and/or statistical forms (including Bayesian minimization, maximum likelihood, and minimum variance methods). The fundamental tenet of data adjustment and data assimilation is the minimization of a user-defined “cost functional” which describes the squared departures between computational results and the observations/experiments of the respective results. In contradistinction to the least-squares tenet underlying data adjustment and/or assimilation, the “BERRU-PM” methodology developed by Cacuci [13] [14] employs the “maximum entropy” (MaxEnt) principle [15] to combine computational and experimental information for obtaining best-estimate predicted mean values for model responses and parameters, together with reduced predicted uncertainties for these best-estimate values, thereby eliminating the need for minimizing the user-chosen “quadratic cost functional representing the weighted errors between measured and computed responses.” BERRU-PM is an acronym for “Best-Estimate Results with Reduced Uncertainties -Predictive Modeling” and is a MaxEnt methodology that incorporates first-order sensitivities of model responses with respect to the model parameters. Such sensitivities are most efficiently computed by using the adjoint sensitivity analysis method for nonlinear systems originally developed by Cacuci [16] [17] . The first-order BERRU-PM methodology was recently extended by Cacuci [18] to include second-order sensitivities. The main differences between DA [5] [6] [7] [8] results and the BERRU-PM [16] [17] [18] results are as follows:

1) DA [5] [6] [7] [8] is formulated conceptually either just in the phase-space of measured responses (“observation-space formulation”) or just in the phase-space of the model’s dependent variables (“state-space formulation”). Hence, DA can calibrate initial conditions as “direct results” but cannot directly calibrate any other model parameters. In contradistinction, the BERRU-PM methodology is formulated conceptually in the most inclusive “joint-phase-space of parameters, computed and measured responses.” Consequently, the BERRU-PM methodology simultaneously calibrates responses and parameters, thus simultaneously providing results for forward and inverse problems.

2) If experiments are perfectly well known, the DA methodology fails fundamentally, but the BERRU-PM methodology does not. The DA methodology also fails fundamentally when the response measurements happen to coincide with the computed value of the response, in which DA yields the trivial result. In contradistinction, the BERRU-PM methodology yields a non-trivial result.

3) The BERRU-PM methodology is significantly more efficient computationally than DA since it only requires the inversion a matrix of size having the dimensions of the number of model responses, whereas DA requires inversion of much larger matrices in the phase-space of dependent variables.

In this work, the methodology presented in [18] will be extended by including the full second-order representation of the computed responses, thereby creating the 2^nd-BERRU-PMD methodology, where the attribute “2^nd” indicates that this methodology incorporates second-order uncertainties (means and covariances) and second-order sensitivities of computed model responses to model parameters, and where the acronym BERRU stands, as before, for “Best-Estimate Results with Reduced Uncertainties;” the last letter (“D”) in the acronym indicates “deterministic,” referring to the deterministic inclusion of the computational model responses. The 2^nd-BERRU-PMD methodology is fundamentally based on the MaxEnt principle.

This work is structured as follows: Section 2 presents the mathematical modeling of the physical system under consideration. Section 3 presents the development of the 2^nd-BERRU-PMD methodology. It is shown that the results produced by the 2^nd-BERRU-PMD methodology include and extend the results produced by the first-order BERRU-PM methodology and the extant data assimilation and/or data adjustment procedures while overcoming the fundamental limitations of these procedures. The concluding discussion presented in Section 4 lays the ground for the subsequent presentation [19] of the alternative framework for developing the “second-order MaxEnt predictive modelling methodology” which probabilistically (as opposed to “deterministically”) incorporates the computed model responses. This alternative methodology [19] will be designated by the acronym 2^nd-BERRU-PMP, where the last letter (“P”) in the acronym indicates “probabilistic” inclusion of the computational model responses.

2. Mathematical Modeling of the Physical System

In general terms, the modeling of a physical system and/or the result of an indirect experimental measurement requires consideration of the following modeling components:

1) A mathematical model comprising linear and/or nonlinear equations that relate the system’s independent variables and parameters to the system’s state (i.e., dependent) variables;

2) Inequality and/or equality constraints that delimit the ranges of the system’s parameters;

3) One or several computational results, customarily referred to as system responses (or objective functions, or indices of performance), which are computed using the mathematical model; and

4) Experimentally measured responses, with their respective nominal (mean) values and uncertainties (variances, covariances, skewness, kurtosis, etc.).

A nonlinear physical system can be generally modeled by means of coupled equations which can be represented in operator form as follows:

$N^m\left[v\left(x\right),\alpha \right]=Q^m\left(x,\alpha \right),x\in {\Omega }_x\left(\alpha \right)$ . (1)

The superscript “m” in Equation (1) indicates “model.” Boundary and/or initial conditions must also be provided if differential operators appear in Equation (1). In operator form, these boundaries and/or initial conditions are represented as follows:

$B^m\left[v\left(x\right);\alpha ;x\right]-B^s\left(x,\alpha \right)=0,x\in \partial {\Omega }_x\left(\alpha \right)$ , (2)

where the vector $B^s\left(x,\alpha \right)$ symbolically indicates “source terms on the boundary.” In this work, matrices will be denoted using capital bold letters while vectors will be denoted using either capital or lower-case bold letters. The symbol “ $\triangleq$ ” will be used to denote “is defined as” or “is by definition equal to.” Transposition will be indicated by a dagger ( $†$ ) superscript. The equalities in this work are considered to hold in the weak (“distributional”) sense. The right-sides of Equations (1) and (2), as well as of other various equations to be derived in this work, may contain “generalized functions/functionals”, particularly Dirac-distributions and derivatives thereof. The quantities which appear in Equations (1) and (2) are defined below.

1) The quantity $\alpha \triangleq {\left({\alpha }_1,\cdots ,{\alpha }_{TP}\right)}^{†}$ is a TP-dimensional vector, having components ${\alpha }_1,\cdots ,{\alpha }_{TP}$ , which denote the model’s imprecisely known parameters. The quantity “TP” denotes the “total number of model parameters.” Without loss of generality, the model parameters $\alpha \triangleq {\left({\alpha }_1,\cdots ,{\alpha }_{TP}\right)}^{†}\in D_{\alpha }\in {ℝ}^{TP}$ can be considered to be real scalars defined over a domain $D_{\alpha }$ , which is included in a TP-dimensional subset of the ${ℝ}^{TP}$ . These model parameters usually stem from processes that are external to the physical system under consideration and their precise value is seldom, if ever, known. The known characteristics of the model parameters may include their nominal (expected/mean) values and, possibly, higher-order moments or cumulants (i.e., variance/covariances, skewness, kurtosis), which are usually determined from experimental data and/or processes external to the physical system under consideration. Occasionally, just the lower and the upper bounds may be known for some model parameters. The components of the TP-dimensional column vector $\alpha \in {ℝ}^{TP}$ are considered to also include imprecisely known geometrical parameters that characterize the physical system’s boundaries in the phase-space of the model’s independent variables. Mathematically, the model parameters can be considered to be quasi-random scalar-valued quantities which follow an unknown multivariate distribution denoted as $p_{\alpha }\left(\alpha \right)$ . The mean values (which in the context of computational modeling are called “nominal” values) of the model parameters will be denoted as ${\alpha }_i^0$ ; the superscript “0” will be used throughout this monograph to denote “nominal values.” These nominal values are formally defined as follows:

${\alpha }_i^0\triangleq {\displaystyle \underset{D_{\alpha }}{\int }{\alpha }_i{p}_{\alpha }\left(\alpha \right)\text{d}\alpha },i=1,\cdots ,TP.$ (3)

The expected values of the measured parameters will be considered to constitute the components of a vector denoted as ${\alpha }^0\triangleq {\left[{\alpha }_1^0,\cdots ,{\alpha }_i^0,\cdots ,{\alpha }_{TP}^0\right]}^{†}$ . The covariance, $\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right)$ , of two model parameters, ${\alpha }_i$ and ${\alpha }_j$ , is defined as follows:

$\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right)\triangleq {\displaystyle \underset{D_{\alpha }}{\int }\left({\alpha }_i-{\alpha }_i^0\right)\left({\alpha }_j-{\alpha }_j^0\right)p_{\alpha }\left(\alpha \right)}\text{d}\alpha ;i,j=1,\cdots ,TP.$ (4)

The covariances $\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right)$ are considered to constitutes the components of the covariance matrix for the model parameters and will be denoted as $C_{\alpha \alpha }\triangleq {\left[\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right)\right]}_{TP\times TP}$ .

2) The generic nonlinear model is considered to comprise TI independent variables which will be denoted as $x_i,i=1,\cdots ,TI$ , and which are considered to be components of a TI-dimensional column vector denoted as $x\triangleq {\left(x_1,\cdots ,x_{TI}\right)}^{†}\in {ℝ}^{TI}$ , where the sub/superscript “TI” denotes the “total number of independent variables.” The vector $x\in {ℝ}^{TI}$ of independent variables is considered to be defined on a phase-space domain which will be denoted as $\Omega \left(\alpha \right)$ and which is defined as follows: $\Omega \left(\alpha \right)\triangleq \left\{-\infty \le {\lambda }_i\left(\alpha \right)\le x_i\le {\omega }_i\left(\alpha \right)\le \infty ;i=1,\cdots ,TI\right\}$ . The lower boundary-point of an independent variable is denoted as ${\lambda }_i\left(\alpha \right)$ and the corresponding upper boundary-point is denoted as ${\omega }_i\left(\alpha \right)$ . The boundary $\Omega \left(\alpha \right)$ is also considered to be imprecisely known since it may depend on both geometrical parameters and material properties. A typical example of boundaries that depend on both geometrical parameters and material properties are the “boundaries facing vacuum” in models based on diffusion theory, where conditions are imposed on the “extrapolated boundary” of the respective spatial domain. The “extrapolated boundary” depends both on the imprecisely known physical dimensions of the problem’s domain and also on the medium’s properties, such as atomic number densities and microscopic transport cross sections. The boundary of $\Omega \left(\alpha \right)$ , denoted as $\partial \Omega \left(\alpha \right)\triangleq \left\{{\lambda }_i\left(\alpha \right)\cup {\omega }_i\left(\alpha \right),i=1,\cdots ,TI\right\}$ , includes the set of the endpoints ${\lambda }_i\left(\alpha \right),{\omega }_i\left(\alpha \right),i=1,\cdots ,TI$ of the respective intervals on which the components of $x$ are defined.

3) The vector $v\left(x\right)\triangleq {\left[v_1\left(x\right),\cdots ,v_{TD}\left(x\right)\right]}^{†}$ is a TD-dimensional column vector, having components $v_i\left(x\right),i=1,\cdots ,TD$ , which represent the model’s dependent variables (also called the model’s “state functions”). The abbreviation “TD” denotes “total number of dependent variables.” Without loss of generality, we can consider that $v\left(x\right)\in E_v$ , where $E_v$ is a normed linear space over the scalar field $F$ of real numbers.

4) The quantity $N^m\left[v\left(x\right);\alpha \right]\triangleq {\left[N_1^m\left(v;\alpha \right),\cdots ,N_{TD}^m\left(v;\alpha \right)\right]}^{†}$ is a TD-dimensional column vector comprising the left-sides of the equations underlying the mathematical/computational model; the superscript “m” indicates “model.” The components $N_i^m\left(v;\alpha \right),i=1,\cdots ,TD$ are operators (including differential, difference, integral, distributions, and/or finite or infinite matrices) acting nonlinearly on the dependent variables $v\left(x\right)$ , the independent variables $x$ and the model parameters $\alpha$ . The mapping $N^m\left(v;\alpha \right)$ is defined on the combined domains of the model’s parameters and state functions, i.e., $N:D\subset E\to E_Q$ , where $D=D_u\oplus D_{\alpha }$ , $D_u\subset E_u$ , $D_{\alpha }\subset E_{\alpha }$ , $E=E_u\oplus E_{\alpha }$ .

5) The vector $Q^m\left(x,\alpha \right)\triangleq {\left[q_1^m\left(x;\alpha \right),\cdots ,q_{TD}^m\left(x;\alpha \right)\right]}^{†}\in E_Q$ is a TD-dimensional column vector which represents inhomogeneous source terms (i.e., the right-sides of the equations underlying the model), which usually depend nonlinearly on the uncertain parameters $\alpha$ . The vector $Q^m\left(x,\alpha \right)$ is defined on a normed linear space denoted as $E_Q$ .

6) In Equation (2), the vector $B^m\left(v;\alpha \right)\triangleq {\left[B_1^m\left(v;\alpha \right),\cdots ,B_{TB}^m\left(v;\alpha \right)\right]}^{†}$ , where the subscript “TB” denotes “total number of boundary conditions,” comprises components $B_i^m\left(v;\alpha \right),i=1,\cdots ,TB$ , which are nonlinear operators in $v\left(x\right)$ and $\alpha$ , which are defined on the boundary $\partial {\Omega }_x\left(\alpha \right)$ of the model’s domain ${\Omega }_x\left(\alpha \right)$ . The vector $B^s\left(x,\alpha \right)\triangleq {\left[B_1^s\left(x;\alpha \right),\cdots ,B_{TB}^s\left(x;\alpha \right)\right]}^{†}$ comprises inhomogeneous boundary sources (indicated by the superscript “s”) which are nonlinear functions of $\alpha$ . The column vector $0$ has TB components, all of which are identically zero.

Solving Equations (1) and (2) at the nominal parameter values, ${\alpha }^0$ , provides the “nominal solution” $v^0\left(x\right)$ , i.e., the vectors $v^0\left(x\right)$ and ${\alpha }^0$ satisfy the following equations:

$N^m\left[v^0\left(x\right);{\alpha }^0\right]=Q^m\left(x,{\alpha }^0\right),x\in {\Omega }_x$ , (5)

$B^m\left[v^0\left(x\right);{\alpha }^0;x\right]-B^s\left(x,{\alpha }^0\right)=0,x\in \partial {\Omega }_x\left({\alpha }^0\right).$ (6)

The results computed using a mathematical model are customarily called “model responses” (or “system responses” or “objective functions” or “indices of performance”). Consider that there are a total number of TR such model responses, each of which can be considered to be a component of the “vector of model responses” $r\triangleq {\left(r_1,\cdots ,r_k,\cdots ,r_{TR}\right)}^{†}$ . Each of these model responses is formally a function (implicit and/or explicit) of the dependent variables and model parameters $\alpha$ , which can be represented formally as follows:

$r_k=r_k\left[v\left(x\right);\alpha \right];\text{}\text{}\text{}k=1,\cdots ,TR.$ (7)

In particular, a measurement of a physical quantity that depends on the model’s state functions and parameters can be considered to be a response denoted as $R_p\left[v\left(x\right);\alpha \right]$ , which is to be evaluated at $x=x^p\left(\alpha \right)$ , where $x^p\left(\alpha \right)\triangleq {\left[x_1^p\left(\alpha \right),\cdots ,x_k^p\left(\alpha \right),\cdots ,x_{TI}^p\left(\alpha \right)\right]}^{†}$ denotes the location in phase-space of the specific “measurement point.” Such a measurement (or measurement-like) response can be represented mathematically as follows:

$\begin{array}{l}{R}_p\left[v\left(x\right);\alpha \right]\\ \triangleq {\displaystyle \underset{{\lambda }_1\left(\alpha \right)}{\overset{{\omega }_1\left(\alpha \right)}{\int }}\cdots {\displaystyle \underset{{\lambda }_{TI}\left(\alpha \right)}{\overset{{\omega }_{TI}\left(\alpha \right)}{\int }}F\left[v\left(x\right);\alpha ;x\right]\delta \left[x_1-x_1^p\left(\alpha \right)\right]}}\times \cdots \times \delta \left[x_{TI}-x_{TI}^p\left(\alpha \right)\right]\text{d}{x}_1\cdots \text{d}{x}_{TI}.\end{array}$ (8)

where the function $F\left[v\left(x\right);\alpha ;x\right]$ comprises the mathematical dependence of the measurement device on the model’s dependent variable(s), and where the quantity $\delta \left[x_i-x_i^p\left(\alpha \right)\right]$ denotes the Dirac-delta functional. The measurement’s location in phase-space, $x^p\left(\alpha \right)$ , may itself be afflicted by measurement (experimental) uncertainties. Hence, the components of $x^p\left(\alpha \right)$ must be included among the components of the vector $\alpha$ of model parameters, even though the quantity $x^p\left(\alpha \right)$ appears only in the definition of the response but does not appear in Equations (1) and (2), which mathematically define the physical model. Thus, the physical “system” is to be understood as being defined to comprise both the system’s computational model and the system’s responses. In most cases, the coordinates $x_k^p\left(\alpha \right)$ , $k=1,\cdots ,TI$ , will simply be independent (albeit uncertain) model parameters included among the components of the vector $\alpha$ , in which case $\partial x_k^p\left(\alpha \right)/\partial {\alpha }_n=1$ , if ${\alpha }_n\equiv x_k^p$ and $\partial x_k^p\left(\alpha \right)/\partial {\alpha }_n=0$ , if ${\alpha }_n\ne x_k^p$ .

The expression on the right-side of Equation (8) can evidently be computed using the model, thus producing the “computed response” value, which will be denoted as $r_k^c\left[v\left(x\right);\alpha ;x\right]$ , and which can be compared to the corresponding “experimentally measured” value; the superscript “c” indicates “computed. The computed response $r_k^c\left[v\left(x\right);\alpha ;x\right]$ depends (implicitly and/or explicitly) on the model’s parameters and dependent variables, which also depend, in turn, on the model’s parameters. Therefore, the uncertainties affecting the model parameters $\alpha$ will “propagate” both directly and indirectly, through the model’s dependent variables, to induce uncertainties in the computed responses, which will therefore be denoted simply as $r_k^c\left(\alpha \right)$ . The nominal value of the response will be denoted $r_k\left({\alpha }^0\right)$ , indicating that this value is obtained by computing the model response using the expected/nominal parameter values ${\alpha }^0\triangleq {\left({\alpha }_1^0,\cdots ,{\alpha }_{TP}^0\right)}^{†}$ .

3. 2^nd-BERRU-PMD: Second Order MaxEnt Predictive Modeling Methodology with Deterministically Included Computed Responses

This Section presents the mathematical and physical considerations leading to the development of the “second-order MaxEnt predictive modeling methodology for obtaining best-estimate results with reduced uncertainties” by deterministically incorporating the computational model. This methodology is designated using the acronym 2^nd-BERRU-PMD, where the last letter (“D”) indicates “deterministic,” referring to the inclusion of the computed model responses. The posterior first-order moments (best-estimate predicted mean values) and posterior second-order moments (best-estimate predicted correlations) produced by the 2^nd-BERRU-PMD methodology are compared with the results produced by the 1^st-BERRU-PM methodology [14] , as well as to the results produced by data assimilation procedures, indicating the specific ways in which the 2^nd-BERRU-PMD overcomes limitations of these latter procedures.

3.1. MaxEnt Representation of Experimental Information for Responses and Parameters

Consider that the total number of experimentally measured responses is TR. The information usually available regarding the distribution of such measured responses comprises the first-order moments (mean values), which will be denoted as $r_i^e$ , $i=1,\cdots ,TR$ , and the second-order moments (variances/covariances), which will be denoted as $\mathrm{cov}{\left(r_i,r_j\right)}_e$ , $i,j=1,\cdots ,TR$ , for the measured responses. The letter “e” will be used either as a superscript or a superscript to indicate experimentally measured quantities. The expected values of the experimentally measured responses will be considered to constitute the components of a vector denoted as $r^e\triangleq {\left(r_1^e,\cdots ,r_{TR}^e\right)}^{†}$ . The covariances of the measured responses are considered to be components of the $TR\times TR$ -dimensional covariance matrix of measured responses, which will be denoted as $C_{rr}^e\triangleq {\left[\mathrm{cov}{\left(r_i,r_j\right)}_e\right]}_{TR\times TR}$ . In principle, it is also possible to obtain correlations between some measured responses and some model parameters. When such correlations between measured responses and measured model parameters are available, they will be denoted as $\text{cor}{\left({\alpha }_i,r_j\right)}_e$ , $i=1,\cdots ,TP;j=1,\cdots ,TR$ , and they can formally be considered to be elements of a rectangular correlation matrix which will be denoted as $C_{\alpha r}^e\triangleq {\left[\text{cor}{\left({\alpha }_i,r_j\right)}_e\right]}_{TP\times TR}$ . As discussed in Section 2, cf. Equations (3) and (4), the model parameters are characterized by the vector of mean values ${\alpha }^0\triangleq {\left[{\alpha }_1^0,\cdots ,{\alpha }_i^0,\cdots ,{\alpha }_{TP}^0\right]}^{†}$ and the covariance matrix $C_{\alpha \alpha }\triangleq {\left[\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right)\right]}_{TP\times TP}$ .

The MaxEnt principle can now be applied, as described in Appendix A, to construct the least informative (and hence, most conservative) distribution using the available experimental information, to obtain the following expression:

$p_e\left(z|z_e,C_e\right)=\frac{{\left(2\pi \right)}^{-\left(TR+TP\right)/2}}{\sqrt{Det\left(C_e\right)}}\exp \left[-\frac{1}{2}{\left(z-z_e\right)}^{†}{C}_e^{-1}\left(z-z_e\right)\right],$ (9)

where $Det\left(C_e\right)$ denotes the determinant of the matrix $C_e$ and where:

$C_e\triangleq \left(\begin{array}{cc}{C}_{rr}^e& C_{r\alpha }^e\\ C_{\alpha r}^e& C_{\alpha \alpha }\end{array}\right);z=\left(\begin{array}{c}r\\ \alpha \end{array}\right);z_e=\left(\begin{array}{c}{r}^e\\ {\alpha }^0\end{array}\right).$ (10)

3.2. Construction of the Joint Posterior MaxEnt Probability Distribution 2^nd-BERRU-PMD

Consistent with the consideration that only mean values and covariances (i.e., first- and second-order distributional moments) are available, the 2^nd-BERRU-PMD methodology considers that only first- and second-order sensitivities of the computed responses with respect to the model parameters are available. In this case, the vector of model responses can be represented by the following second-order Taylor-series written in “vector/matrix” form:

$\begin{array}{l}r=r^c\left({\alpha }^0\right)+\left[S\left({\alpha }^0\right)\right]\left(\alpha -{\alpha }^0\right)+p\left(\alpha ;{\alpha }^0\right);\\ p\left(\alpha ;{\alpha }^0\right)\triangleq {\left[p_1\left(\alpha ;{\alpha }^0\right),\cdots ,p_{TR}\left(\alpha ;{\alpha }^0\right)\right]}^{†};\\ p_k\left(\alpha ;{\alpha }^0\right)\triangleq \frac{1}{2}{\left(\alpha -{\alpha }^0\right)}^{†}\left[Q_k\left({\alpha }^0\right)\right]\left(\alpha -{\alpha }^0\right),\end{array}$ (11)

where the components of the vector $r^c\left({\alpha }^0\right)\triangleq {\left[r_1^c\left({\alpha }^0\right),\cdots ,r_{TR}^c\left({\alpha }^0\right)\right]}^{†}$ represent the values of the model responses computed at the nominal parameter values ${\alpha }^0$ , while the other vectors and matrices contain first-order and second-order sensitivities of responses with respect to the model parameters, evaluated at the nominal parameter values ${\alpha }^0$ , and are defined below.

$\begin{array}{l}S\left(\alpha \right)\triangleq \left(\begin{array}{ccc}{s}_{11}& \cdots & s_{1,TR}\\ ⋮& \ddots & ⋮\\ s_{TR,1}& \cdots & s_{TR,TR}\end{array}\right)\triangleq \left(\begin{array}{ccc}\frac{\partial r_1\left(\alpha \right)}{\partial {\alpha }_1}& \cdots & \frac{\partial r_1\left(\alpha \right)}{\partial {\alpha }_{TP}}\\ ⋮& \ddots & ⋮\\ \frac{\partial r_{TR}\left(\alpha \right)}{\partial {\alpha }_1}& \cdots & \frac{\partial r_{TR}\left(\alpha \right)}{\partial {\alpha }_{TP}}\end{array}\right);\\ Q_k\left(\alpha \right)\triangleq \left(\begin{array}{ccc}{q}_{11}^{\left(k\right)}& \cdots & q_{1,TP}^{\left(k\right)}\\ ⋮& \ddots & ⋮\\ q_{TP,1}^{\left(k\right)}& \cdots & q_{TP,TP}^{\left(k\right)}\end{array}\right)\triangleq \left(\begin{array}{ccc}\frac{{\partial }^2{r}_k^c\left(\alpha \right)}{\partial {\alpha }_1\partial {\alpha }_1}& \cdots & \frac{{\partial }^2{r}_k^c\left(\alpha \right)}{\partial {\alpha }_1\partial {\alpha }_{TP}}\\ ⋮& \ddots & ⋮\\ \frac{{\partial }^2{r}_k^c\left(\alpha \right)}{\partial {\alpha }_{TP}\partial {\alpha }_1}& \cdots & \frac{{\partial }^2{r}_k^c\left(\alpha \right)}{\partial {\alpha }_{TP}\partial {\alpha }_{TP}}\end{array}\right).\end{array}$ (12)

The merging of the experimental information with the computational information using the mathematical model is accomplished by seeking a mathematical representation of a “true” probability distribution $p\left(z\right)$ , $z\triangleq {\left(r,\alpha \right)}^{†}$ , which incorporates all of the second-order information provided in the foregoing about the model parameters, measured and computed responses. This “true” probability distribution must satisfy the following characteristics: 1) $p\left(z\right)$ is normalized to unity; 2) $p\left(z\right)$ is approximated as closely as possible by the experimental distribution $p_e\left(z|z^e,C_e\right)$ and 3) $p\left(z\right)$ satisfies the following integral forms of the deterministic relationships among responses and parameters represented mathematically by Equation (11):

$\begin{array}{l}{\displaystyle \underset{Dz}{\int }{g}_k\left(z\right)p\left(z\right)\text{d}}z=0;\\ g_k\left(z\right)\triangleq r_k-r_k^c\left({\alpha }^0\right)-{\displaystyle \underset{i=1}{\overset{TP}{\sum }}{s}_{ki}\delta {\alpha }_i}-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{TP}{\sum }}{\displaystyle \underset{j=1}{\overset{TP}{\sum }}{q}_{ij}^{\left(k\right)}}}\delta {\alpha }_i\delta {\alpha }_j=0;k=1,\cdots ,TR.\end{array}$ (13)

The “true” normalized distribution $p\left(z\right)$ which minimizes the discrepancy between it and the distribution $p_e\left(z|z^e,C_e\right)$ while satisfying the constraints expressed by Equation (13) is obtained by applying the steps generally outlined in Appendix B, which yields the following expressions:

$p\left(z\right)=\frac{p_e\left(z|z^e,C_e\right)}{Z}\exp \left\{-{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k{g}_k\left(z\right)}\right\};$ (14)

$Z\triangleq {\displaystyle \underset{Dz}{\int }{p}_e\left(z|z^e,C_e\right)\exp \left\{-{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k{g}_k\left(z\right)}\right\}\text{d}}z.$ (15)

In Equations (14) and (15), the quantities ${\theta }_1,\cdots ,{\theta }_k,\cdots ,{\theta }_{TR}$ are Lagrange multipliers which remain to be determined. In statistical mechanics, the normalization integral Z defined in Equation (15) is called the partition function (or sum over states) and carries all of the information available about the possible states of the system. The quantity $\left(-\log Z\right)$ is called the “free energy” and plays an important role in determining the expressions of the Lagrange multipliers and hence of the probability distribution $p\left(z\right)$ . In the present case, it follows from Equations (15) and (13) that:

$\frac{\partial \left[-\ln Z\right]}{\partial {\theta }_k}=\frac{1}{Z}{\displaystyle \underset{Dz}{\int }{g}_k\left(z\right)p_e\left(z|z^e,C_e\right)\exp \left\{-{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k{g}_k\left(z\right)}\right\}\text{d}z}=0,k=1,\cdots ,TR.$ (16)

On the other hand, alternative expressions for the derivatives $\partial \left[-\ln Z\right]/\partial {\theta }_k$ of the free energy with respect to the Lagrange multipliers can be obtained by evaluating the integral which defines the partition function Z in Equation (15). Equating these alternative expressions to zero, in view of Equation (16), will provide equations which are to be solved to obtain the explicit expressions of the Lagrange multipliers.

Integrals such as represented by Equation (15) can be quantified, to any desired order of accuracy—even exactly—using the well-known saddle-point (Laplace expansion) method, see e.g., [20] . Examining the exponent of the function entering the normalization integral Z defined in Equation (15) readily indicates that if the term involving the second-order sensitivities of the response with respect to the parameters is retained in the exponent, the expression of the saddle-point will involve the vector of Lagrange multipliers, $\theta \triangleq {\left({\theta }_1,\cdots ,{\theta }_k,\cdots ,{\theta }_{TR}\right)}^{†}$ , in a nonlinear equation which cannot be solved analytically to obtain a closed-form expression for these Lagrange multipliers. Furthermore, solving numerically this nonlinear equation in $\theta$ would require the inversion of the (very large) parameter-covariance matrix $C_{\alpha \alpha }$ , which is impractical. On the other hand, as will be shown below, a closed-form expression for the saddle-point and, consequently, for the Lagrange multipliers can be obtained if the terms containing the second-order response sensitivities are considered by expanding the respective exponential term in a Taylor-series in powers of $\theta \triangleq {\left({\theta }_1,\cdots ,{\theta }_k,\cdots ,{\theta }_{TR}\right)}^{†}$ or, equivalently, around the nominal parameter values. Expanding the second-order terms transforms the normalization integral Z defined in Equation (15) into the following form:

$Z\triangleq {\displaystyle \underset{Dz}{\int }f\left(r,\alpha \right){\text{e}}^{-h\left(r,\alpha ;\theta \right)}\text{d}r\text{d}\alpha }$ , (17)

where $\theta \triangleq {\left({\theta }_1,\cdots ,{\theta }_k,\cdots ,{\theta }_{TR}\right)}^{†}$ and where:

$\begin{array}{c}h\left(r,\alpha ;\theta \right)\triangleq \frac{1}{2}{\left(\begin{array}{c}r-r^e\\ \alpha -{\alpha }^0\end{array}\right)}^{†}{\left(\begin{array}{cc}{C}_{rr}^e& C_{r\alpha }^e\\ C_{\alpha r}^e& C_{\alpha \alpha }\end{array}\right)}^{-1}\left(\begin{array}{c}r-r^e\\ \alpha -{\alpha }^0\end{array}\right)\\ +{\theta }^{†}\left[-r+r^c\left({\alpha }^0\right)+S\left(\alpha -{\alpha }^0\right)\right],\end{array}$ (18)

$f\left(r,\alpha \right)\triangleq 1-\frac{1}{2}{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k}{\left(\alpha -{\alpha }^0\right)}^{†}\left[Q_k\left({\alpha }^0\right)\right]\left(\alpha -{\alpha }^0\right)+hots,$ (19)

where the acronym “hots” stands for “higher-order terms in the expansion of the exponential $\exp \left\{-{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k{g}_k\left(z\right)}\right\}$ .” Neglecting these higher-order terms, it follows from Equations (17)-(19) that

$Z=I_1-\frac{1}{2}{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k}{J}_k$ , (20)

where

$I_1\triangleq {\displaystyle \underset{Dz}{\int }{\text{e}}^{-h\left(r,\alpha ;\theta \right)}\text{d}r\text{d}\alpha }$ , (21)

$J_k\triangleq {\displaystyle \underset{Dz}{\int }{\left(\alpha -{\alpha }^0\right)}^{†}\left[Q_k\left({\alpha }^0\right)\right]\left(\alpha -{\alpha }^0\right){\text{e}}^{-h\left(r,\alpha ;\theta \right)}\text{d}r\text{d}\alpha }$ . (22)

The saddle-point, denoted as $z^{sp}\triangleq {\left(r^{sp},{\alpha }^{sp}\right)}^{†}$ , of the function $h\left(r,\alpha ;\theta \right)$ is defined as the point where the following partial gradients vanish:

${\nabla }_{r}h\left(r,\alpha ;\theta \right)=0,{\nabla }_{\alpha }h\left(r,\alpha ;\theta \right)=0,\text{at}z=z^{sp}.$ (23)

The conditions represented by Equation (23) yield the following equation:

${\left(\begin{array}{cc}{C}_{rr}^e& C_{r\alpha }^e\\ C_{\alpha r}^e& C_{\alpha \alpha }\end{array}\right)}^{-1}\left(\begin{array}{c}{r}^{sp}-r^e\\ {\alpha }^{sp}-{\alpha }^0\end{array}\right)=\left(\begin{array}{c}\theta \\ -S^{†}\theta \end{array}\right).$ (24)

Solving Equation (24) yields the following expressions for the components of the saddle point $z^{sp}\triangleq {\left(r^{sp},{\alpha }^{sp}\right)}^{†}$ of the function $h\left(r,\alpha ;\theta \right)$ :

$r^{sp}=r^e+\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right)\theta$ , (25)

${\alpha }^{sp}={\alpha }^0+\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right)\theta$ . (26)

Using the saddle point method (which is exact for the integral $I_1$ ) or “completing the square in the exponent of $I_1$ ” yields the following exact expression for the integral $I_1$ :

$I_1=\left(2\pi \right){\text{}}^{\left(TR+TP\right)/2}{\left|C_e\right|}^{1/2}\exp \left[-{\theta }^{†}d+\frac{1}{2}{\left(\begin{array}{c}\theta \\ -S^{†}\theta \end{array}\right)}^{†}{C}_e\left(\begin{array}{c}\theta \\ -S^{†}\theta \end{array}\right)\right]$ , (27)

where $\left|C_e\right|$ denotes the determinant of the matrix $C_e$ and where $d$ denotes the “vector of deviations between the nominal values of computed and measured responses” and is defined below:

$d\triangleq r^c\left({\alpha }^0\right)-r^e;\text{}d\triangleq {\left(d_1,\cdots ,d_{TR}\right)}^{†};d_1\triangleq r_i^c\left({\alpha }^0\right)-r_i^e;i=1,\cdots ,TR.$ (28)

The expression of the response in terms of the computational model, i.e., Equation (11), can be used to obtain the following representation for $d$ in terms of the computed responses:

$\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}d\triangleq r-r^e-S\left(\alpha -{\alpha }^0\right)-p\left(\alpha ;{\alpha }^0\right)$ . (29)

Thus, the “vector of deviations” $d$ can be considered to be a (quasi-) random variable. The expressions of the “vector of mean values of deviations” $\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\bar{d}\triangleq {\left[{\bar{d}}_1,\cdots ,{\bar{d}}_{TR}\right]}^{†}$ and the covariance matrix $D\triangleq {\left[D_{ij}\right]}_{TR\times TR}$ of the vector of deviations $d$ are derived below:

$\text{}\text{}\text{}\text{}\text{}\text{}\text{}\text{}\bar{d}\triangleq {\displaystyle \underset{Dz}{\int }\left[r-r^e-S\left(\alpha -{\alpha }^0\right)-p\left(\alpha ;{\alpha }^0\right)\right]p_e\left(z|z^e,C_e\right)\text{d}z}={\left[{\bar{d}}_1,\cdots ,{\bar{d}}_{TR}\right]}^{†}$ , (30)

${\bar{d}}_k\triangleq {\displaystyle \underset{i=1}{\overset{TP}{\sum }}{\displaystyle \underset{j=1}{\overset{TP}{\sum }}{q}_{ij}^{\left(k\right)}\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right)}};k=1,\cdots ,TR$ ; (31)

$\begin{array}{c}D\triangleq {\displaystyle \underset{Dz}{\int }\left(d{d}^{†}\right)p_e\left(z|z^e,C_e\right)\text{d}z}\\ ={\displaystyle \underset{Dz}{\int }{p}_e\left(z|z^e,C_e\right)\left[r-r^e-S\left(\alpha -{\alpha }^0\right)-p\left(\alpha ;{\alpha }^0\right)\right]}\\ \times {\left[r-r^e-S\left(\alpha -{\alpha }^0\right)-p\left(\alpha ;{\alpha }^0\right)\right]}^{†}\text{d}z\\ \cong C_{rr}^e+S{C}_{\alpha \alpha }{S}^{†}-C_{r\alpha }^e{S}^{†}-S{C}_{\alpha r}^e.\end{array}$ (32)

The third- and fourth-order correlations were neglected in the final expression obtained for the covariance matrix $D\triangleq {\left[D_{ij}\right]}_{TR\times TR}$ in Equation (32). Notably, the result obtained in Equation (30) indicates that if the second-order sensitivities of the response with respect to the model parameters are neglected, then ${\bar{d}}_k=0,k=1,\cdots ,TR$ , in which case $d$ would be a random vector with zero-mean and covariance $D$ . Therefore, neglecting the second-order sensitivities would imply that the computed model responses would be consistent with the measurements only if $\bar{d}=0$ , i.e., $r^c\left({\alpha }^0\right)\cong r^e$ . Although such a situation would not be impossible to conceive, it would probably not occur often in practice. These considerations further underscore the significant impact which the second-order response sensitivities have on the correct interpretation of actual, physical systems.

Applying the saddle-point method to the expression provided in Equation (22) yields the following expression:

$J_k=I_1{\left({\alpha }^{sp}-{\alpha }^0\right)}^{†}\left[Q_k\left({\alpha }^0\right)\right]\left({\alpha }^{sp}-{\alpha }^0\right)\left[1+O\left(C_{\alpha \alpha }\right)\right],$ (33)

where the quantity $O\left(C_{\alpha \alpha }\right)$ denotes terms of the order of the components of the parameter covariance matrix $C_{\alpha \alpha }$ , and can therefore be neglected by comparison to the leading term of order unity.

Replacing the expression obtained in Equations (27) and (33) into Equation (20) and neglecting the higher-order terms yields the following expression for the partition function:

$Z=I_1\left\{1-\frac{1}{2}{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k}{\left({\alpha }^{sp}-{\alpha }^0\right)}^{†}\left[Q_k\left({\alpha }^0\right)\right]\left({\alpha }^{sp}-{\alpha }^0\right)\right\}$ . (34)

Taking the logarithm of Equation (34) yields the following expression for the “free energy”:

$-\ln Z=\frac{TR+TP}{2}\ln \left(2\pi \right)+\frac{1}{2}\ln \text{}\left|C_m\right|-{\theta }^{†}d+h_1\left(\theta \right)+h_2\left(\theta \right)$ , (35)

$h_1\left(\theta \right)\triangleq \frac{1}{2}\left[{\theta }^{†}{C}_{rr}^e\theta -{\theta }^{†}\left(C_{r\alpha }^e{S}^{†}+S{C}_{\alpha r}^e\right)\theta +{\theta }^{†}S{C}_{\alpha \alpha }{S}^{†}\theta \right]$ , (36)

$\begin{array}{c}\text{}{h}_2\left(\theta \right)\triangleq \ln \left\{1-\frac{1}{2}{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k}{\left({\alpha }^{sp}-{\alpha }^0\right)}^{†}\left[Q_k\left({\alpha }^0\right)\right]\left({\alpha }^{sp}-{\alpha }^0\right)\right\}\\ =\ln \left\{1-\frac{1}{2}{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k\left[{\theta }^{†}{M}^{\left(k\right)}\theta \right]}\right\},\end{array}$ (37)

$M^{\left(k\right)}\triangleq {\left[M_{ij}^{\left(k\right)}\right]}_{TR\times TR}\triangleq {\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right)}^{†}\left[Q_k\left({\alpha }^0\right)\right]\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right)={\left[M^{\left(k\right)}\right]}^{†}$ . (38)

The expression for $h_2\left(\theta \right)$ presented in Equation (37) has been obtained by using Equation (26). As indicated in Equations (36) and (37), respectively, the quantity $h_1\left(\theta \right)$ contains only first-order response sensitivities (with respect to the model’s parameters) and second-order bilinear terms in the components of $\theta$ , while the quantity $h_2\left(\theta \right)$ comprises all second-order response sensitivities and terms higher than second-order in the components of $\theta$ .

In view of Equations (35)-(37), the derivatives of the free energy with respect to the Lagrange multipliers ${\theta }_k$ are provided by the expression below:

$\frac{\partial \left[-\ln Z\right]}{\partial {\theta }_k}=-d_k+\frac{\partial h_1\left(\theta \right)}{\partial {\theta }_k}+\frac{\partial h_2\left(\theta \right)}{\partial {\theta }_k},k=1,\cdots ,TR$ . (39)

It follows from Equations (39) and (16) that the vector of Lagrange multipliers must be the solution of the following system of coupled nonlinear equations:

$w_k\left(\theta \right)\triangleq \frac{\partial h_1\left(\theta \right)}{\partial {\theta }_k}+\frac{\partial h_2\left(\theta \right)}{\partial {\theta }_k}-d_k=0,k=1,\cdots ,TR,$ (40)

where:

$\frac{\partial h_1\left(\theta \right)}{\partial {\theta }_k}={\displaystyle \underset{j=1}{\overset{TR}{\sum }}{D}_{kj}}{\theta }_j;{\left[D_{kj}\right]}_{TR\times TR}\triangleq D\triangleq C_{rr}^e+S{C}_{\alpha \alpha }{S}^{†}-C_{r\alpha }^e{S}^{†}-S{C}_{\alpha r}^e;$ (41)

$\text{}\frac{\partial h_2\left(\theta \right)}{\partial {\theta }_k}=-\frac{1}{2}\frac{\frac{\partial }{\partial {\theta }_k}{\displaystyle \underset{i=1}{\overset{TR}{\sum }}{\theta }_i\left[{\theta }^{†}{M}^{\left(i\right)}\theta \right]}}{1-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{TR}{\sum }}{\theta }_i\left[{\theta }^{†}{M}^{\left(i\right)}\theta \right]}}\cong -\frac{1}{2}{\theta }^{†}{M}^{\left(k\right)}\theta -{\displaystyle \underset{i=1}{\overset{TR}{\sum }}{\theta }_i}{\displaystyle \underset{j=1}{\overset{TR}{\sum }}{M}_{kj}^{\left(i\right)}{\theta }_j}\triangleq u_k\left(\theta \right).$ (42)

The approximation $1-\frac{1}{2}{\displaystyle \underset{j=1}{\overset{TR}{\sum }}{\theta }_j\left[{\theta }^{†}{M}^{\left(j\right)}\theta \right]}\cong 1$ has been used to obtain the approximate expression on the right-side of Equation (42), which thus comprises just the bilinear terms in $\theta$ , neglecting the high-order terms involving the components of $\theta$ , which are expected to be comparatively small.

Introducing the vectors $w\left(\theta \right)\triangleq {\left[w_1\left(\theta \right),\cdots ,w_k\left(\theta \right),\cdots ,w_{TR}\left(\theta \right)\right]}^{†}$ and $u\left(\theta \right)\triangleq {\left[u_1\left(\theta \right),\cdots ,u_k\left(\theta \right),\cdots ,u_{TR}\left(\theta \right)\right]}^{†}$ enables the compact writing, in matrix-vector form, of Equation (40), as follows:

$w\left(\theta \right)\triangleq -d+D\theta +u\left(\theta \right)=0$ . (43)

If only the first-order sensitivities of the response with respect to the model parameters are considered (and thus the second- and higher-order sensitivities are neglected or are unavailable), the term $u\left(\theta \right)$ is neglected and the remaining first-order term in Equation (43) yields the following first-order solution, denoted as ${\theta }^{\left(1\right)}$ , for the vector of Lagrange multipliers:

${\theta }^{\left(1\right)}=D^{-1}d.$ (44)

Notably, the matrix $D$ has dimensions $TR\times TR$ . Since $TR\ll TP$ in practical situations, it follows that the computation of the coordinates of the saddle-point $z^{sp}\triangleq {\left(r^{sp},{\alpha }^{sp}\right)}^{†}$ involves the inversion of just the smallest-possible matrix (i.e., $D$ ) that arises in the “predictive modeling” methodology. Replacing the result obtained for ${\theta }^{\left(1\right)}$ in Equations (25) and (26) provides the corresponding expressions for the coordinates $r^{sp}$ and ${\alpha }^{sp}$ of the saddle-point $z^{sp}\triangleq {\left(r^{sp},{\alpha }^{sp}\right)}^{†}$ .

Using the results obtained for ${\theta }^{\left(1\right)}$ in Equation leads to the expressions already presented in the 1^st-BERRU-PM methodology [14] for the saddle point $z^{sp}\triangleq {\left(r^{sp},{\alpha }^{sp}\right)}^{†}$ and, subsequently for the best-estimate (calibrated) mean values for the parameters and responses, as well as for the best-estimate covariance matrices for the corresponding best-estimate parameters and responses. These expressions will not be reproduced here, since they will also be obtained as particular cases of the 2^nd-BERRU-PM methodology to be presented in the remainder of this Section.

The goal of predictive modeling is to reduce, as much possible, the discrepancies between the measured response values and the corresponding values predicted by the computational model But if the nominal values of the measured responses happen—by nature, design or accident—to coincide with the corresponding computed values, the vector of deviations vanishes, i.e., $d=0$ , in which case Equations (44), (25), (26) and (14) indicate that:

$d=0;\Rightarrow {\theta }^{\left(1\right)}=0;\Rightarrow r^{sp}=r^e,{\alpha }^{sp}={\alpha }^0;\Rightarrow p\left(z\right)=p_e\left(z|z^e,C_e\right).$ (45)

The results shown in Equation (45) indicate that if only first-order sensitivities are used and if the nominal values of the measured responses happen to coincide with the corresponding computed values, the model would have no impact on the initial (measured) distribution of responses and parameters, so there would be no calibration of either the predicted responses or parameters. Evidently this conclusion reveals a fundamental limitation of using just the first-order sensitivities for predictive modeling (and data assimilation and/or data adjustment methods). But this conclusion also provides the fundamental motivation for using second-order sensitivities, in addition to the natural expectation that using the second-order sensitivities will produce more accurate results, in general.

It is also important to note that if the response measurements were perfect and if they were uncorrelated to the model parameters, i.e., if $C_{rr}^e=0$ and $C_{r\alpha }^e=C_{\alpha r}^e=0$ , the matrix $D=S{C}_{\alpha \alpha }{S}^{†}$ would nevertheless be non-singular since the matrix $S{C}_{\alpha \alpha }{S}^{†}$ is nonsingular. This, the solution shown in Equation (44) would still exist; this is in contradistinction with the extant data assimilation schemes, which fail if the measurements are perfect. Notably, all extant data assimilation schemes consider that the response measurements and model parameters are uncorrelated.

If the second-order sensitivities, which are included in the term $u\left(\theta \right)$ , are retained in Equation (40), this equation becomes nonlinear in the components of $\theta \triangleq {\left({\theta }_1,\cdots ,{\theta }_k,\cdots ,{\theta }_{TR}\right)}^{†}$ and therefore requires an iterative method for solving. The structure of $u\left(\theta \right)$ indicates that the simplest (and probably the most efficient) method for solving Equation (43) would be Newton’s method, as can be surmised by the similarity between this equation and a simple quadratic algebraic equation. In this vein, consider the application of Newton’s method to the quadratic equation, as described below:

$\begin{array}{l}q\left(x\right)\triangleq a{x}^2+bx+d=0;q^{\prime }\left(x\right)=2ax+b;\\ x^{\left(n+1\right)}=x^{\left(n\right)}-{\left[q^{\prime }\left(x^n\right)\right]}^{-1}q\left(x^n\right);\text{}\text{}n=0,1,\cdots .\end{array}$ (46)

If $d\ne 0$ , the usual initial-guess for starting the Newton iteration is $x^{\left(0\right)}=-d/b$ , which yields the first-iterate $x^{\left(1\right)}=-\left[1+ad/\left(b^2-2ad\right)\right]d/b$ . The correction term included in this first-iterate solution is of first order in the product ( $ad$ ), so the subsequent iterates become increasingly more accurate when this product is small; even the first-iterate becomes exact in the limit as $a\to 0$ . On the other hand, if $d=0$ , the initial guess $x^{\left(0\right)}=-d/b=0$ would yield solely the trivial solution $x=0$ . Therefore, if $d\cong 0$ , the initial-guess for starting the Newton iteration would be $x^{\left(0\right)}=-b/a$ , which would yield the non-zero first-iterate solution $x^{\left(1\right)}=-b/a+d/b$ , and so on.

Applying Newton’s method to solve Equation (43) yields the following iterative equation for finding its solution(s) $\theta \triangleq {\left({\theta }_1,\cdots ,{\theta }_k,\cdots ,{\theta }_{TR}\right)}^{†}$ :

${\theta }^{\left(n+1\right)}={\theta }^{\left(n\right)}-{\left[A\left({\theta }^{\left(n\right)}\right)\right]}^{-1}w\left({\theta }^{\left(n\right)}\right);\text{}\text{}n=0,1,\cdots ,$ (47)

where the Jacobian matrix $A\left(\theta \right)$ is defined as usual, namely:

$\begin{array}{c}A\left(\theta \right)\triangleq \frac{\partial \left(w_1,\cdots ,w_k,\cdots ,w_{TR}\right)}{\partial \left({\theta }_1,\cdots ,{\theta }_n,\cdots ,{\theta }_{TR}\right)}\triangleq \left[\begin{array}{ccc}\partial w_1/\partial {\theta }_1& ·& \partial w_1/\partial {\theta }_{TR}\\ ·& ·& ·\\ \partial w_{TR}/\partial {\theta }_1& ·& \partial w_{TR}/\partial {\theta }_{TR}\end{array}\right]\\ \triangleq \left[\begin{array}{ccc}{A}_{1,1}\left(\theta \right)& ·& A_{1,TR}\left(\theta \right)\\ ·& ·& ·\\ A_{TR,1}\left(\theta \right)& ·& A_{TR,TR}\left(\theta \right)\end{array}\right].\end{array}$ (48)

Using Equations (42) and (43) yields the following expression for the components/elements $A_{k,n}\left(\theta \right)$ of the Jacobian matrix $A\left(\theta \right)$ :

$\begin{array}{c}{A}_{kn}\left(\theta \right)=\frac{\partial w_k}{\partial {\theta }_n}=\frac{\partial }{\partial {\theta }_n}\left[{\displaystyle \underset{j=1}{\overset{TR}{\sum }}{D}_{kj}}{\theta }_j-\frac{1}{2}{\theta }^{†}{M}^{\left(k\right)}\theta -{\displaystyle \underset{i=1}{\overset{TR}{\sum }}{\theta }_i}{\displaystyle \underset{j=1}{\overset{TR}{\sum }}{M}_{kj}^{\left(i\right)}{\theta }_j}\right]\\ =D_{kn}-{\displaystyle \underset{j=1}{\overset{TR}{\sum }}{M}_{nj}^{\left(k\right)}{\theta }_j}-{\displaystyle \underset{j=1}{\overset{TR}{\sum }}{M}_{kj}^{\left(n\right)}{\theta }_j}-{\displaystyle \underset{i=1}{\overset{TR}{\sum }}{\theta }_i{M}_{kn}^{\left(i\right)}}=D_{kn}-V_{kn}\left(\theta \right),\end{array}$ (49)

where:

$V_{kn}\left(\theta \right)\triangleq {\displaystyle \underset{j=1}{\overset{TR}{\sum }}{\theta }_j\left[M_{nj}^{\left(k\right)}+M_{kj}^{\left(n\right)}+M_{kn}^{\left(j\right)}\right]},V\left(\theta \right)\triangleq {\left[V_{kn}\left(\theta \right)\right]}_{TR\times TR}$ . (50)

As indicated in Equation (49), if the initial guess ${\theta }^{\left(0\right)}=0$ is used, the Jacobian matrix would take on the value $A\left(\theta =0\right)=D$ . Consequently, the resulting solution of Equation (47) would be obtained in a single iteration and would be the same as the first-order solution obtained in Equation (44), i.e., ${\theta }^{\left(1\right)}=D^{-1}d$ . Similarly, if the second-order response sensitivities to the model parameters are neglected, Equation (50) indicates that $A\left(M^{\left(k\right)}=0\right)=D$ .

To obtain the nontrivial solution of Equation (47) when $d$ is very small or actually vanishes ( $d=0$ ), a non-zero initial guess ( ${\theta }^{\left(0\right)}\ne 0$ ) must be used for starting the Newton iterations, as was discussed in the foregoing by analogy with the solution of a quadratic algebraic equation using Newton’s method. The nontrivial solution of Equation (47) which would be valid even when $d=0$ will be denoted as ${\theta }^{\left(h,n+1\right)}$ , where the superscript “h” indicates that the expression of ${\theta }^{\left(h,n+1\right)}$ will include higher-order terms, which will involve the second-order sensitivities of computed responses with respect to model parameters. The corresponding Newton iteration indicated by Equation (47) will thus take on the following form for finding the nontrivial root of Equation (43) when $d$ is very small or actually vanishes:

${\theta }^{\left(h,n+1\right)}={\theta }^{\left(h,n\right)}-{\left[D-V\left({\theta }^{\left(h,n\right)}\right)\right]}^{-1}w\left({\theta }^{\left(h,n\right)}\right);\text{}\text{}n=0,1,\cdots .$ (51)

As has been discussed in the foregoing in conjunction with the algebraic quadratic equation, it would be advantageous to start the Newton iteration shown in Equation (51) with an initial guess ${\theta }^{\left(h,0\right)}$ that would have components analogous to the initial guess $x^{\left(0\right)}=-b/a$ for the quadratic equations. It would be further advantageous if the computation of the initial guess ${\theta }^{\left(h,0\right)}$ would not require inversion of matrices. An expression for ${\theta }^{\left(h,0\right)}$ which satisfies these considerations can be derived by setting $d=0$ in Equation (43) and by considering just the diagonal terms of the matrices involved in this equation, thus obtaining:

${\theta }^{\left(h,0\right)}\triangleq {\left({\theta }_1^{\left(h,0\right)},\cdots ,{\theta }_j^{\left(h,0\right)},\cdots ,{\theta }_{TR}^{\left(h,0\right)}\right)}^{†};{\theta }_j^{\left(h,0\right)}=2D_{jj}/3M_{jj}^{\left(j\right)};j=1,\cdots ,TR.$ (52)

Using the initial guess provided in Equation (52) in Equation (50) yields the following starting expressions for components of the Jacobian matrix $A\left(\theta \right)$ :

$A_{kn}\left({\theta }^{\left(h,0\right)}\right)=D_{kn}-{\displaystyle \underset{j=1}{\overset{TR}{\sum }}{\theta }_j^{\left(h,0\right)}\left[M_{nj}^{\left(k\right)}+M_{kj}^{\left(n\right)}+M_{kn}^{\left(j\right)}\right]}.$ (53)

Furthermore, using the initial guess provided in Equation (52) in Equation (43) yields the following starting expressions for components $w_k\left(\theta \right)$ of the vector $w\left(\theta \right)\triangleq {\left[w_1\left(\theta \right),\cdots ,w_k\left(\theta \right),\cdots ,w_{TR}\left(\theta \right)\right]}^{†}$ :

$w_k\left({\theta }^{\left(h,0\right)}\right)=u_k\left({\theta }^{\left(h,0\right)}\right)+{\displaystyle \underset{j=1}{\overset{TR}{\sum }}{D}_{kj}}{\theta }_j^{\left(h,0\right)}-d_k,k=1,\cdots ,TR;$ (54)

$\begin{array}{c}{u}_k\left({\theta }^{\left(h,0\right)}\right)=-\frac{1}{2}{\displaystyle \underset{i=1}{\overset{TR}{\sum }}{\displaystyle \underset{j=1}{\overset{TR}{\sum }}{\theta }_i^{\left(h,0\right)}{M}_{ij}^{\left(k\right)}{\theta }_j^{\left(h,0\right)}}}-{\displaystyle \underset{i=1}{\overset{TR}{\sum }}{\theta }_i^{\left(h,0\right)}}{\displaystyle \underset{j=1}{\overset{TR}{\sum }}{M}_{kj}^{\left(i\right)}{\theta }_j^{\left(h,0\right)}}\\ \cong -\frac{4}{9}{\displaystyle \underset{i=1}{\overset{TR}{\sum }}{\displaystyle \underset{j=1}{\overset{TR}{\sum }}\left(\frac{D_{ii}}{M_{ii}^{\left(i\right)}}\right)\left(\frac{D_{jj}}{M_{jj}^{\left(j\right)}}\right)\left[\frac{1}{2}{M}_{ij}^{\left(k\right)}+M_{kj}^{\left(i\right)}\right]}},k=1,\cdots ,TR.\end{array}$ (55)

The impact of the second-order sensitivities through the components $M_{ij}^{\left(k\right)}$ of the matrices $M^{\left(k\right)}$ is evident in Equations (53)-(55). The first-iterate ${\theta }^{\left(h,1\right)}$ of the solution of Equation (51) has the following expression:

$\begin{array}{l}{\theta }^{\left(h,1\right)}={\theta }^{\left(h,0\right)}-{\left(D-V^{\left(h,0\right)}\right)}^{-1}\left(u^{\left(h,0\right)}+D{\theta }^{\left(h,0\right)}-d\right);\\ V^{\left(h,0\right)}\triangleq V\left({\theta }^{\left(h,0\right)}\right);u^{\left(h,0\right)}\triangleq u\left({\theta }^{\left(h,0\right)}\right).\end{array}$ (56)

It is evident from the expressions provided in Equations (53)-(55) that even if $d=0$ (i.e., if the vector of deviations between the nominally measured and computed values of the responses vanishes), the starting-values of the components $w_k\left({\theta }^{\left(h,0\right)}\right)\ne 0$ are non-zero because of the impact of the second-order response sensitivities. In addition, the initial-Jacobian matrix $A\left({\theta }^{\left(h,0\right)}\right)$ is nonsingular (also because of the impact of the second-order response sensitivities), so the first-iterate ${\theta }^{\left(h,1\right)}$ has a finite non-zero value. All subsequent iterates will also have non-zero values. Thus, using the second-order response sensitivities overcomes the limitation of the first-order predictive modeling methodology 1^st-BERRU-PM, as well as the limitation of extant data adjustment and data assimilation procedures, which all fail when the measured response values coincide with the computed response values.

Replacing the solution of Equation (51) into Equations (25) and (26) yields the following expressions for the coordinates of the saddle point $z^{sp}\triangleq {\left(r^{sp},{\alpha }^{sp}\right)}^{†}$ of the function $h\left(r,\alpha ;\theta \right)$ :

$r^{sp}=r^e+\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right){\theta }^{\left(h,n+1\right)}$ , (57)

${\alpha }^{sp}={\alpha }^0+\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right){\theta }^{\left(h,n+1\right)}$ . (58)

The best-estimate (optimal) mean values (first-order moments) and correlations (second-order moments) of the posterior distribution of the combined computational model with measured responses and parameters are obtained by using the saddle-point method to evaluate the integrals that define these first- and second-order moments. These best-estimate values will be denoted using the superscript “bed” to indicate that these “best-estimate” values arise from deterministically incorporated/assimilated computed responses.

Thus, the predicted “best-estimated values with deterministically incorporated computed responses” for the predicted mean values of responses and parameters, respectively, have the following expressions:

$r^{bed}\triangleq {\displaystyle \underset{Dz}{\int }rp\left(z\right)\text{d}z}=r^{sp}=r^e+\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right){\theta }^{\left(h,n+1\right)}$ ; (59)

${\alpha }^{bed}\triangleq {\displaystyle \underset{Dz}{\int }\alpha p\left(z\right)\text{d}z}={\alpha }^{sp}={\alpha }^0+\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right){\theta }^{\left(h,n+1\right)}$ . (60)

The expression of the “best-estimate with deterministically incorporated computed responses” covariance matrix for the responses, denoted as $C_{rr}^{bed}$ , is obtained from its definition in conjunction with Equation (59), which yields:

$\begin{array}{l}{C}_{rr}^{bed}\triangleq {\displaystyle \underset{Dz}{\int }\left(r-r^{bed}\right){\left(r-r^{bed}\right)}^{†}p\left(z\right)\text{d}z}\\ ={\displaystyle \underset{Dz}{\int }\left[r-r^e-\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right){\theta }^{\left(h,n+1\right)}\right]{\left[r-r^e-\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right){\theta }^{\left(h,n+1\right)}\right]}^{†}p\left(z\right)\text{d}z}.\end{array}$ (61)

The expression of the “best-estimate with deterministically incorporated computed responses” covariance matrix for the parameters, denoted as $C_{\alpha \alpha }^{bed}$ , is obtained from its definition in conjunction with Equation (60), which yields:

$\begin{array}{l}{C}_{\alpha \alpha }^{bed}\triangleq {\displaystyle \underset{Dz}{\int }\left(\alpha -{\alpha }^{bed}\right){\left(\alpha -{\alpha }^{bed}\right)}^{†}p\left(z\right)\text{d}z}\\ ={\displaystyle \underset{Dz}{\int }\left[\alpha -{\alpha }^0-\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right){\theta }^{\left(h,n+1\right)}\right]{\left[\alpha -{\alpha }^0-\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right){\theta }^{\left(h,n+1\right)}\right]}^{†}p\left(z\right)\text{d}z}.\end{array}$ (62)

The expression of the “best-estimate with deterministically incorporated computed responses” correlation matrix among parameters and responses, denoted as $C_{\alpha r}^{bed}$ , is obtained by using its definition in conjunction with Equations (59) and (60), which yields:

$\begin{array}{l}{C}_{\alpha r}^{bed}\triangleq {\displaystyle \underset{Dz}{\int }\left(\alpha -{\alpha }^{bed}\right){\left(r-r^{bed}\right)}^{†}p\left(z\right)\text{d}z}\\ ={\displaystyle \underset{Dz}{\int }\left[\alpha -{\alpha }^0-\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right){\theta }^{\left(h,n+1\right)}\right]{\left[r-r^e-\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right){\theta }^{\left(h,n+1\right)}\right]}^{†}p\left(z\right)\text{d}}z\text{}.\end{array}$ (63)

The impact of the second-order response sensitivities on the mean values and correlation/covariance matrices can be assessed already by using the first iterate solution, ${\theta }^{\left(h,1\right)}$ , in Equations (59)-(63). Using the first-iterate, ${\theta }^{\left(h,1\right)}$ , in Equations (59) and (60), respectively, yields the following expressions for the “first-iterate best-estimate” values for the responses and parameters, respectively:

$r^{\left(1\right)bed}=r^e+\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right)\left[{\theta }^{\left(h,0\right)}-{\left(D-V^{\left(h,0\right)}\right)}^{-1}\left(u^{\left(h,0\right)}+D{\theta }^{\left(h,0\right)}-d\right)\right]\text{}\text{}\text{}\text{}\text{}\text{}$ ; (64)

${\alpha }^{\left(1\right)bed}={\alpha }^0+\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right)\left[{\theta }^{\left(h,0\right)}-{\left(D-V^{\left(h,0\right)}\right)}^{-1}\left(u^{\left(h,0\right)}+D{\theta }^{\left(h,0\right)}-d\right)\right]$ . (65)

As the expressions obtained in Equations (64) and (65) indicate, combining experimental and computational information modifies the correlations among the parameters through couplings introduced by the sensitivities of the participating responses. The expressions obtained in Equations (64) and (65) possess the following new features by comparisons to their counterparts provided by extant data assimilation and data adjustment procedures:

1) If the second-order sensitivities are neglected, then the expressions in (64) and (65) take on the following simplified forms:

$\text{no}{\text{2}}^{\text{nd}}\text{-order}\text{sensitivities}:r^{\left(1\right)bed}=r^e+\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right)D^{-1}d$ ; (66)

$\text{no}{\text{2}}^{\text{nd}}\text{-order}\text{sensitivities}:{\alpha }^{\left(1\right)bed}={\alpha }^0+\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right)D^{-1}d$ . (67)

As expected, the expressions in Equations (66) and (67) are identical to the expressions provided by the 1^st-BERRU-PM [14] methodology and therefore possess the same strengths (easy to implement) and weaknesses (trivial results, no calibration if $d=0$ ) as the results produced by the 1^st-BERRU-PM methodology.

2) If the vector of deviations between the nominally measured and computed values of the responses vanishes, i.e., if $d=0$ , the expressions provided in Equations (64) and (65) take on the following simplified forms:

$\begin{array}{l}\text{if}d=0:\\ r^{\left(1\right)bed}\left(d=0\right)=r^e+\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right)\left[{\theta }^{\left(h,0\right)}-{\left(D-V^{\left(h,0\right)}\right)}^{-1}\left(u^{\left(h,0\right)}+D{\theta }^{\left(h,0\right)}\right)\right]\ne 0\end{array}$ (68)

$\begin{array}{l}\text{if}d=0:\\ {\alpha }^{\left(1\right)bed}\left(d=0\right)={\alpha }^0+\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right)\left[{\theta }^{\left(h,0\right)}-{\left(D-V^{\left(h,0\right)}\right)}^{-1}\left(u^{\left(h,0\right)}+D{\theta }^{\left(h,0\right)}\right)\right]\ne 0\end{array}$ (69)

As evidenced by comparing Equations (68) and (69) with Equations (66) and (67), respectively, the contributions of the second-order sensitivities provide the fundamental difference when $d=0$ between non-calibration [i.e., $r^{\left(1\right)bed}\left(d=0\right)=r^e$ and ${\alpha }^{\left(1\right)bed}\left(d=0\right)={\alpha }^0$ ] in the case for data adjustment, data assimilation and 1^st-BERRU-PM methodologies, cf. Equations (66) and (67), respectively, versus the calibration $r^{\left(1\right)bed}\left(d=0\right)\ne r^e$ and ${\alpha }^{\left(1\right)bed}\left(d=0\right)\ne {\alpha }^0$ provided by the 2^nd-BERRU-PMD expressions obtained in Equations (68) and (69).

3) The components of the matrix $V^{\left(h,0\right)}$ are usually smaller (since they contain higher-order quantities) than the corresponding components of the matrix $D$ , so the contributions provided by the term $V^{\left(h,0\right)}$ can be considered to be “corrections” to the components of the matrix $D$ . If the contributions of the matrix $V^{\left(h,0\right)}$ are neglected by comparison to the corresponding components of the matrix $D$ , the expressions provided in Equations (64) and (65) reduce to the following expressions:

$\text{if}{V}^{\left(h,0\right)}=0:r^{\left(1\right)bed}\left(V^{\left(h,0\right)}=0\right)=r^e+\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right)D^{-1}\left(d-u^{\left(h,0\right)}\right)\text{}\text{}\text{}\text{}\text{}\text{}$ ; (70)

$\text{if}{V}^{\left(h,0\right)}=0:{\alpha }^{\left(1\right)bed}\left(V^{\left(h,0\right)}=0\right)={\alpha }^0+\left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right)D^{-1}\left(d-u^{\left(h,0\right)}\right)$ . (71)

The above expressions underscore, through the term $u^{\left(h,0\right)}$ , the fundamental importance of including second-order sensitivities, since they provide the basic contribution to maintaining the effects of calibration for both responses and parameters, even when $d=0$ , i.e., when the values of the measured responses coincide with the corresponding values of the computed responses. This essential contribution stemming from the second-order responses sensitivities is clearly highlighted by comparing the expressions obtained in Equations (70) and (71) with the corresponding first-order (1^st-BERRU-PM) expressions provided in Equations (66) and (67), respectively.

Using the first-iterate, ${\theta }^{\left(h,1\right)}$ , in Equation (61) yields the following expression for the “first-iterate best-estimate” covariance matrix $C_{rr}^{\left(1\right)bed}$ for the responses:

$\begin{array}{l}{C}_{rr}^{\left(1\right)bed}={\displaystyle \underset{Dz}{\int }\left\{r-r^e-\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right)\left[{\theta }^{\left(h,0\right)}-{\left(D-V^{\left(h,0\right)}\right)}^{-1}\left(u^{\left(h,0\right)}+D{\theta }^{\left(h,0\right)}-d\right)\right]\right\}}\text{}\\ \times {\left\{r-r^e-\left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right)\left[{\theta }^{\left(h,0\right)}-{\left(D-V^{\left(h,0\right)}\right)}^{-1}\left(u^{\left(h,0\right)}+D{\theta }^{\left(h,0\right)}-d\right)\right]\right\}}^{†}p\left(z\right)\text{d}z\\ =\langle \left[r-r^e+X_r{\theta }^{\left(h,0\right)}-Y_{r}d+Y_r{u}^{\left(h,0\right)}\right]{\left[r-r^e+X_r{\theta }^{\left(h,0\right)}-Y_{r}d+Y_r{u}^{\left(h,0\right)}\right]}^{†}\rangle .\end{array}$ (72)

The following definitions were used in Equation (72):

$\langle •\rangle \triangleq {\displaystyle \underset{Dz}{\int }(•)p\left(z\right)\text{d}z}\text{},$ (73)

$\begin{array}{c}{C}_{rd}^{\left(1\right)}\triangleq \langle \left(r-r^e\right)d^{†}\rangle =\langle \left(r-r^e\right){\left[r-r^e-S\left(\alpha -{\alpha }^0\right)-p\left(\alpha ;{\alpha }^0\right)\right]}^{†}\rangle \\ \cong C_{rr}^e-C_{r\alpha }^e{S}^{†};\end{array}$ (74)

$X_r\triangleq C_{rd}^{\left(1\right)}\left[{\left(D-V^{\left(h,0\right)}\right)}^{-1}D-I\right]$ ; (75)

$Y_r\triangleq C_{rd}^{\left(1\right)}{\left(D-V^{\left(h,0\right)}\right)}^{-1}$ . (76)

Expanding the integrand in the second equality in Equation (72) yields the following expression for the “first-iterate best-estimate” response covariance matrix $C_{rr}^{\left(1\right)bed}$ :

$\begin{array}{c}{C}_{rr}^{\left(1\right)bed}=C_{rr}^e-C_{rd}^{\left(1\right)}{Y}_r^{†}+X_r{\theta }^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_r^{†}-X_r{\theta }^{\left(h,0\right)}{\bar{d}}^{†}{Y}_r^{†}+X_r{\theta }^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}\\ -Y_r{\left(C_{rd}^{\left(1\right)}\right)}^{†}-Y_r\bar{d}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_r^{†}+Y_{r}D{Y}_r^{†}-Y_r\bar{d}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}\\ +Y_r{u}^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_r^{†}-Y_r{u}^{\left(h,0\right)}{\bar{d}}^{†}{Y}_r^{†}+Y_r{u}^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}.\end{array}$ (77)

Using the first-iterate, ${\theta }^{\left(h,1\right)}$ , in Equation (62) and evaluating the resulting integrals yields the following expression for the “first-iterate best-estimate” covariance matrix $C_{\alpha \alpha }^{\left(1\right)be}$ for the model parameters:

$\begin{array}{c}{C}_{\alpha \alpha }^{\left(1\right)bed}=\langle \left[\alpha -{\alpha }^0+X_{\alpha }{\theta }^{\left(h,0\right)}-Y_{\alpha }d+Y_{\alpha }{u}^{\left(h,0\right)}\right]{\left[\alpha -{\alpha }^0+X_{\alpha }{\theta }^{\left(h,0\right)}-Y_{\alpha }d+Y_{\alpha }{u}^{\left(h,0\right)}\right]}^{†}\rangle \\ =C_{\alpha \alpha }^e-C_{\alpha d}^{\left(1\right)}{Y}_{\alpha }^{†}+X_{\alpha }{\theta }^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_{\alpha }^{†}-X_{\alpha }{\theta }^{\left(h,0\right)}{\bar{d}}^{†}{Y}_{\alpha }^{†}+X_{\alpha }{\theta }^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_{\alpha }^{†}\\ -Y_{\alpha }{\left(C_{\alpha d}^{\left(1\right)}\right)}^{†}-Y_{\alpha }\bar{d}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_{\alpha }^{†}+Y_{\alpha }D{Y}_{\alpha }^{†}-Y_{\alpha }\bar{d}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_{\alpha }^{†}\\ +Y_{\alpha }{u}^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_{\alpha }^{†}-Y_{\alpha }{u}^{\left(h,0\right)}{\bar{d}}^{†}{Y}_{\alpha }^{†}+Y_{\alpha }{u}^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_{\alpha }^{†}.\end{array}$ (78)

The following definitions were used in Equation (78):

$\begin{array}{c}{C}_{\alpha d}^{\left(1\right)}\triangleq \langle \left(\alpha -{\alpha }^0\right)d^{†}\rangle =\langle \left(\alpha -{\alpha }^0\right){\left[r-r^e-S\left(\alpha -{\alpha }^0\right)-p\left(\alpha ;{\alpha }^0\right)\right]}^{†}\rangle \\ \cong C_{\alpha r}^e-C_{\alpha \alpha }^e{S}^{†};\end{array}$ (79)

$X_{\alpha }\triangleq C_{\alpha d}^{\left(1\right)}\left[{\left(D-V^{\left(h,0\right)}\right)}^{-1}D-I\right]$ ; (80)

$Y_{\alpha }\triangleq C_{\alpha d}^{\left(1\right)}{\left(D-V^{\left(h,0\right)}\right)}^{-1}$ . (81)

Using the first-iterate, ${\theta }^{\left(h,1\right)}$ , in Equation (63) and evaluating the resulting integrals yields the following expression for the “first-iterate best-estimate” correlation matrix $C_{\alpha r}^{\left(1\right)bed}$ for the parameters and responses:

$\begin{array}{c}{C}_{\alpha r}^{\left(1\right)bed}=\langle \left[\alpha -{\alpha }^0+X_{\alpha }{\theta }^{\left(h,0\right)}-Y_{\alpha }d+Y_{\alpha }{u}^{\left(h,0\right)}\right]{\left[r-r^e+X_r{\theta }^{\left(h,0\right)}-Y_{r}d+Y_r{u}^{\left(h,0\right)}\right]}^{†}\rangle \\ =C_{\alpha r}^e-C_{\alpha d}^{\left(1\right)}{Y}_r^{†}+X_{\alpha }{\theta }^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_r^{†}-X_{\alpha }{\theta }^{\left(h,0\right)}{\bar{d}}^{†}{Y}_r^{†}+X_{\alpha }{\theta }^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}\\ -Y_{\alpha }{\left(C_{rd}^{\left(1\right)}\right)}^{†}-Y_{\alpha }\bar{d}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_r^{†}+Y_{\alpha }D{Y}_r^{†}-Y_{\alpha }\bar{d}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}\\ +Y_{\alpha }{u}^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_r^{†}-Y_{\alpha }{u}^{\left(h,0\right)}{\bar{d}}^{†}{Y}_r^{†}+Y_{\alpha }{u}^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}.\end{array}$ (82)

As the expressions obtained in Equations (77), (78) and (82) indicate, combining experimental and computational information modifies the correlations among the parameters through couplings introduced by the sensitivities of the participating responses.

The expression obtained in Equations (77), (78) and (82) for $C_{rr}^{\left(1\right)bed}$ , $C_{\alpha \alpha }^{\left(1\right)bed}$ and $C_{\alpha r}^{\left(1\right)bed}$ , respectively, will now be analyzed below by considering the same simplifying conditions as were considered for the best-estimate response and parameter mean-values $r^{\left(1\right)bed}$ and ${\alpha }^{\left(1\right)bed}$ , respectively, as follows:

1) If all of the second-order sensitivities are neglected, then the following simplifications occur: $X_r=X_{\alpha }=0$ , $Y_r=C_{rd}^{\left(1\right)}{D}^{-1}$ and $Y_{\alpha }=C_{\alpha d}^{\left(1\right)}{D}^{-1}$ . Therefore, the expressions in Equations (77), (78) and (82) take on the following simplified forms:

$\text{no}{\text{2}}^{\text{nd}}\text{-order}\text{sensitivities}:\text{}\text{}\text{}\text{}\text{}\text{}{C}_{rr}^{\left(1\right)bed}=C_{rr}^e-C_{rd}^{\left(1\right)}{D}^{-1}{C}_{rd}^{\left(1\right)}{}^{†}$ . (83)

$\text{no}{\text{2}}^{\text{nd}}\text{-order}\text{sensitivities}:\text{}\text{}\text{}\text{}\text{}\text{}{C}_{\alpha \alpha }^{\left(1\right)bed}=C_{\alpha \alpha }^e-C_{\alpha d}^{\left(1\right)}{D}^{-1}{C}_{\alpha d}^{\left(1\right)}{}^{†}$ . (84)

$\text{no}{\text{2}}^{\text{nd}}\text{-order}\text{sensitivities}:\text{}\text{}\text{}\text{}\text{}\text{}{C}_{\alpha r}^{\left(1\right)bed}=C_{\alpha r}^e-C_{\alpha d}^{\left(1\right)}{D}^{-1}{C}_{rd}^{\left(1\right)}{}^{†}$ . (85)

As expected, the expressions obtained in Equation (83)-(85) are identical to the expressions provided by the first-order predictive modeling methodology 1^st-BERRU-PM [14] and therefore possess the same strengths (easy to implement) and weaknesses (i.e., trivial results and no calibration if $d=0$ ) as 1^st-BERRU-PM.

2) If the vector of deviations between the nominally measured and computed values of the responses vanishes, i.e., if $d=0$ , the expressions obtained in Equations (77) for $C_{rr}^{\left(1\right)bed}$ , $C_{\alpha \alpha }^{\left(1\right)bed}$ and $C_{\alpha r}^{\left(1\right)bed}$ , respectively, will take on the following simplified forms:

$\begin{array}{l}\text{if}d=0:C_{rr}^{\left(1\right)bed}=C_{rr}^e-C_{rd}^{\left(1\right)}{Y}_r^{†}+X_r{\theta }^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_r^{†}+X_r{\theta }^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}\\ -Y_r{\left(C_{rd}^{\left(1\right)}\right)}^{†}+Y_{r}D{Y}_r^{†}+Y_r{u}^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_r^{†}+Y_r{u}^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†};\end{array}$ (86)

$\begin{array}{l}\text{if}d=0:C_{\alpha \alpha }^{\left(1\right)bed}=C_{\alpha \alpha }^e-C_{\alpha d}^{\left(1\right)}{Y}_{\alpha }^{†}+X_{\alpha }{\theta }^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_{\alpha }^{†}+X_{\alpha }{\theta }^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_{\alpha }^{†}\\ -Y_{\alpha }{\left(C_{\alpha d}^{\left(1\right)}\right)}^{†}+Y_{\alpha }D{Y}_{\alpha }^{†}+Y_{\alpha }{u}^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_{\alpha }^{†}+Y_{\alpha }{u}^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_{\alpha }^{†};\end{array}$ (87)

$\begin{array}{l}\text{if}d=0:C_{\alpha r}^{\left(1\right)bed}=C_{\alpha r}^e-C_{\alpha d}^{\left(1\right)}{Y}_r^{†}+X_{\alpha }{\theta }^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_r^{†}+X_{\alpha }{\theta }^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}\\ -Y_{\alpha }{\left(C_{rd}^{\left(1\right)}\right)}^{†}+Y_{\alpha }D{Y}_r^{†}+Y_{\alpha }{u}^{\left(h,0\right)}{\left({\theta }^{\left(h,0\right)}\right)}^{†}{X}_r^{†}+Y_{\alpha }{u}^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}.\end{array}$ (88)

If the contributions of the second-order sensitivities are neglected in the expressions obtained in Equations (86)-(88), these expressions will reduce to the first-order expressions shown in Equations (83)-(85).

3) If the contributions of the matrix $V^{\left(h,0\right)}$ are neglected by comparison to the corresponding components of the matrix $D$ , then the following simplifications occur: $X_r=X_{\alpha }=0$ , $Y_r=C_{rd}^{\left(1\right)}{D}^{-1}$ and $Y_{\alpha }=C_{\alpha d}^{\left(1\right)}{D}^{-1}$ . Consequently, the expressions provided in Equations (77), (78) and (82) will reduce to the following forms:

$\begin{array}{l}\text{if}{V}^{\left(h,0\right)}=0:C_{rr}^{\left(1\right)bed}\left(V^{\left(h,0\right)}=0\right)=C_{rr}^e-C_{rd}^{\left(1\right)}{D}^{-1}{\left(C_{rd}^{\left(1\right)}\right)}^{†}\\ -Y_r\bar{d}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}-Y_r{u}^{\left(h,0\right)}{\bar{d}}^{†}{Y}_r^{†}+Y_r{u}^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†};\end{array}$ (89)

$\begin{array}{l}\text{if}{V}^{\left(h,0\right)}=0:C_{\alpha \alpha }^{\left(1\right)bed}\left(V^{\left(h,0\right)}=0\right)=C_{\alpha \alpha }^e-C_{\alpha d}^{\left(1\right)}{D}^{-1}{\left(C_{\alpha d}^{\left(1\right)}\right)}^{†}\\ -Y_{\alpha }\bar{d}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_{\alpha }^{†}-Y_{\alpha }{u}^{\left(h,0\right)}{\bar{d}}^{†}{Y}_{\alpha }^{†}+Y_{\alpha }{u}^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_{\alpha }^{†};\end{array}$ (90)

$\begin{array}{l}\text{if}{V}^{\left(h,0\right)}=0:C_{\alpha r}^{\left(1\right)bed}\left(V^{\left(h,0\right)}=0\right)=C_{\alpha r}^e-C_{\alpha d}^{\left(1\right)}{D}^{-1}{\left(C_{rd}^{\left(1\right)}\right)}^{†}\\ -\text{}{Y}_{\alpha }\bar{d}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}-Y_{\alpha }{u}^{\left(h,0\right)}{\bar{d}}^{†}{Y}_r^{†}+Y_{\alpha }{u}^{\left(h,0\right)}{\left(u^{\left(h,0\right)}\right)}^{†}{Y}_r^{†}.\end{array}$ (91)

The above expressions underscore, through the terms involving the vector $u^{\left(h,0\right)}$ , the fundamental importance of including second-order sensitivities, which provide non-zero contributions even when $\bar{d}=0$ .

It is important to emphasize that the posterior distribution of responses and parameters is non-Gaussian since the “best estimate” triple correlations among the best-estimate responses and parameters are nonzero even if the experimentally measured responses and parameters were normally distributed initially. This fact can readily be shown by considering Equations (64) and (65), which can be written as follows:

$\begin{array}{l}{r}^{\left(1\right)bed}=r^e+r^{(0)};\\ r^{\left(0\right)}\triangleq \left(C_{rr}^e-C_{r\alpha }^e{S}^{†}\right)\left[{\theta }^{\left(h,0\right)}-{\left(D-V^{\left(h,0\right)}\right)}^{-1}\left(u^{\left(h,0\right)}+D{\theta }^{\left(h,0\right)}-d\right)\right];\end{array}$ (92)

$\begin{array}{l}{\alpha }^{\left(1\right)bed}={\alpha }^0+{\beta }^{\left(0\right)};\\ {\beta }^{\left(0\right)}\triangleq \left(C_{\alpha r}^e-C_{\alpha \alpha }{S}^{†}\right)\left[{\theta }^{\left(h,0\right)}-{\left(D-V^{\left(h,0\right)}\right)}^{-1}\left(u^{\left(h,0\right)}+D{\theta }^{\left(h,0\right)}-d\right)\right].\end{array}$ (93)

Using Equations (92) and (93) in conjunction with the definition of triple correlations for multivariate random quantities yields the following “lowest-order estimate” for the best-estimate triple correlations among responses and parameters, which denoted below as $T_{rr,ijk}^{\left(1\right)bed}$ and $T_{\alpha \alpha ,ijk}^{\left(1\right)bed}$ , respectively:

$\begin{array}{l}{T}_{rr,ijk}^{\left(1\right)bed}\triangleq \langle \left(r_i-r_i^e-r_i^{\left(0\right)}\right)\left(r_j-r_j^e-r_j^{\left(0\right)}\right)\left(r_k-r_k^e-r_k^{\left(0\right)}\right)\rangle \\ =E\left[\left(r_i-r_i^e\right)\left(r_j-r_j^e\right)\left(r_k-r_k^e\right)\right]-r_i^{\left(0\right)}\mathrm{cov}\left(r_j,r_k\right)-r_j^{\left(0\right)}\mathrm{cov}\left(r_i,r_k\right)\\ -r_k^{\left(0\right)}\mathrm{cov}\left(r_i,r_j\right)-r_i^{\left(0\right)}{r}_j^{\left(0\right)}{r}_k^{\left(0\right)};i,j,k=1,\cdots ,TR;\end{array}$ (94)

$\begin{array}{l}{T}_{\alpha \alpha ,ijk}^{\left(1\right)bed}\triangleq \langle \left({\alpha }_i-{\alpha }_i^0-{\beta }_i^{\left(0\right)}\right)\left({\alpha }_j-{\alpha }_j^0-{\beta }_j^{\left(0\right)}\right)\left({\alpha }_k-{\alpha }_k^e-{\beta }_k^{\left(0\right)}\right)\rangle \\ =E\left[\left({\alpha }_i-{\alpha }_i^0\right)\left({\alpha }_j-{\alpha }_j^0\right)\left({\alpha }_j-{\alpha }_j^0\right)\right]-{\beta }_i^{\left(0\right)}\mathrm{cov}\left({\alpha }_j,{\alpha }_k\right)\\ -{\beta }_j^{\left(0\right)}\mathrm{cov}\left({\alpha }_i,{\alpha }_k\right)-{\beta }_k^{\left(0\right)}\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right)-{\beta }_i^{\left(0\right)}{\beta }_j^{\left(0\right)}{\beta }_k^{\left(0\right)};i,j,k=1,\cdots ,TP.\end{array}$ (95)

If the responses are normally distributed initially, then the initial triple correlations are all zero, i.e., $E\left[\left({\alpha }_i-{\alpha }_i^0\right)\left({\alpha }_j-{\alpha }_j^0\right)\left({\alpha }_j-{\alpha }_j^0\right)\right]$ for all $i,j,k=1,\cdots ,TR$ . Nevertheless, as indicated in Equation (94), the “best-estimate response triple correlations” $T_{rr,ijk}^{\left(1\right)bed}\ne 0$ are non-zero even in such a case, unless there is “no calibration,” i.e., if $r_i^{\left(0\right)}=r_j^{\left(0\right)}=r_k^{\left(0\right)}=0$ for all $i,j,k=1,\cdots ,TR$ . Similarly, if the parameters are normally distributed initially, then the initial triple correlations are all zero, i.e., $E\left[\left(r_i-r_i^e\right)\left(r_j-r_j^e\right)\left(r_k-r_k^e\right)\right]=0$ for all $i,j,k=1,\cdots ,TP$ . Nevertheless, as indicated in Equation (95), the “best-estimate parameter triple correlations” $T_{\alpha \alpha ,ijk}^{\left(1\right)bed}\ne 0$ are non-zero even in such a case, unless there is “no calibration,” i.e., if ${\alpha }_i^{\left(0\right)}={\alpha }_j^{\left(0\right)}={\alpha }_k^{\left(0\right)}=0$ for all $i,j,k=1,\cdots ,TP$ .

4. Discussion and Conclusions

This work has presented the 2^nd-BERRU-PMD methodology, where: the attribute “2^nd” indicates that this methodology incorporates second-order uncertainties (means and covariances) and second-order sensitivities of computed model responses to model parameters; the acronym BERRU stands for “best-estimate results with reduced uncertainties;” and the letter “D” indicates “deterministic,” referring to the deterministic inclusion of the computational model responses. The 2^nd-BERRU-PMD methodology is fundamentally based on the maximum entropy (MaxEnt) principle. This principle is in contradistinction to the fundamental principle that underlies the extant data assimilation and/or adjustment procedures which minimize in a least-square sense a subjective user-defined functional which is meant to represent the discrepancies between measured and computed model responses. It has shown that the 2^nd-BERRU-PMD methodology includes all of the previous results produced by the first-order BERRU-PM methodology, generalizing and extending current data assimilation and/or data adjustment procedures while overcoming the fundamental limitations of these procedures. Since the framework of the 2^nd-BERRU-PMD methodology comprises the combined phase-space of parameters and responses, the expressions obtained for the predicted best-estimate responses and calibrated parameters can be used for solving both forward/direct and inverse problems.

The accompanying work [19] will present the alternative framework for developing the “second-order MaxEnt predictive modelling methodology” by incorporating the computed model responses probabilistically (as opposed to “deterministically”). This alternative methodology [19] will be designated by the acronym 2^nd-BERRU-PMP, where the last letter (“P”) in the acronym indicates “probabilistic” inclusion of the computational model responses. It will be shown that although both the 2^nd-BERRU-PMD and the 2^nd-BERRU-PMP frameworks predict analogous 2^nd-order expressions for the best-estimate posterior first- and second-order moments, the respective expressions for the predicted responses, calibrated predicted parameters and their predicted uncertainties (covariances) are not identical to each other. The additional considerations that make the corresponding expressions predicted by the 2^nd-BERRU-PMD and the 2^nd-BERRU-PMP become identical will also be analyzed in the accompanying Part II [19] .

Appendix A: Construction of the Second-Order Maximum Entropy Distribution for Responses and Parameters

When an unknown distribution $p\left(\alpha ,r\right)$ , defined on a domain $D\triangleq D_{\alpha }\cup D_r$ (where $D_{\alpha }$ denotes the domain of definition of the parameters and $D_r$ denotes the domain of definition of the responses) needs to be reconstructed from a finite number of its known moments, the principle of maximum entropy (MaxEnt) originally formulated by Jaynes [15] provides the optimal compatibility with the available information, while simultaneously ensuring minimal spurious information content. In particular, when only the first-order and second-order moments of the joint distribution of model parameters and responses $p\left(\alpha ,r\right)$ are known, the corresponding second-order MaxEnt distribution $p\left(\alpha ,r\right)$ is constructed by following the procedure outlined below.

1) Known means and covariances for responses:

$r_i^0\triangleq {\displaystyle \underset{D}{\int }{r}_{i}p\left(\alpha ,r\right)\text{d}}r\text{d}\alpha ,i=1,\cdots ,TR;$ (96)

$\mathrm{cov}\left(r_i,r_j\right)\triangleq {\displaystyle \underset{D_e}{\int }\left(r_i-r_i^0\right)\left(r_j-r_j^0\right)p\left(\alpha ,r\right)\text{d}}r\text{d}\alpha ;i,j=1,\cdots ,TR.$ (97)

2) Known means and covariances for parameters:

${\alpha }_j^0\triangleq {\displaystyle \underset{D_e}{\int }{\alpha }_{j}p\left(\alpha ,r\right)\text{d}}r\text{d}\alpha ,j=1,\cdots ,TP;$ (98)

$\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right)\triangleq {\displaystyle \underset{D_e}{\int }\left({\alpha }_i-{\alpha }_i^0\right)\left({\alpha }_j-{\alpha }_j^0\right)p\left(\alpha ,r\right)\text{d}}r\text{d}\alpha ;i,j=1,\cdots ,TP;$ (99)

$\mathrm{cov}\left({\alpha }_i,r_j\right)\triangleq {\displaystyle \underset{D_e}{\int }\left({\alpha }_i-{\alpha }_i^0\right)\left(r_j-r_j^0\right)p\left(\alpha ,r\right)\text{d}}r\text{d}\alpha ;i=1,\cdots ,TP;j=1,\cdots ,TR.$ (100)

According to the MaxEnt principle, the probability density $p\left(\alpha ,r\right)$ would satisfy the “available information” provided in Equations (96)-(100), without implying any spurious information or hidden assumptions, if:

1) $p\left(\alpha ,r\right)$ maximizes the Shannon [21] information entropy, S, as defined below:

$S=-{\displaystyle \underset{D}{\int }p\left(\alpha ,r\right)\ln \left[p\left(\alpha ,r\right)\right]\text{d}\alpha \text{d}r}$ , (101)

2) $p\left(\alpha ,r\right)$ satisfies the “moments constraints” defined by Equations (96)-(100);

3) $p\left(\alpha ,r\right)$ satisfies the normalization condition:

${\displaystyle \underset{D}{\int }p\left(\alpha ,r\right)\text{d}\alpha \text{d}r}=1.$ (102)

The MaxEnt distribution $p\left(\alpha ,r\right)$ is obtained as the solution of the constrained variational problem $\partial H\left(p\right)/\partial p=0$ , where the entropy (Lagrangian functional) $H\left(p\right)$ is defined as follows:

$\begin{array}{c}H\left(p\right)=-{\displaystyle \underset{D}{\int }p\left(\alpha ,r\right)\ln \left[p\left(\alpha ,r\right)\right]\text{d}\alpha \text{d}r}-{\lambda }_0\left[{\displaystyle \underset{D}{\int }p\left(\alpha ,r\right)\text{d}\alpha \text{d}r}-1\right]\\ -{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\lambda }_k^{\left(1\right)}\left[{\displaystyle \underset{D}{\int }{r}_{k}p\left(\alpha ,r\right)\text{d}\alpha \text{d}r-r_k^0}\right]}-{\displaystyle \underset{i=1}{\overset{TP}{\sum }}{\lambda }_i^{\left(2\right)}\left[{\displaystyle \underset{D}{\int }{\alpha }_{i}p\left(\alpha ,r\right)\text{d}\alpha \text{d}r}-{\alpha }_i^0\right]}\end{array}$

$\begin{array}{c}-\frac{1}{2}{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\displaystyle \underset{\mathcal{l}=1}{\overset{TR}{\sum }}{\lambda }_{k\mathcal{l}}^{\left(11\right)}}\left[{\displaystyle \underset{D}{\int }{r}_k{r}_{\mathcal{l}}p\left(\alpha ,r\right)\text{d}\alpha \text{d}r-cov\left(r_k,r_{\mathcal{l}}\right)-r_k^0{r}_{\mathcal{l}}^0}\right]}\\ -{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\displaystyle \underset{i=1}{\overset{TP}{\sum }}{\lambda }_{ki}^{\left(12\right)}}\left[{\displaystyle \underset{D}{\int }{r}_k{\alpha }_{i}p\left(\alpha ,r\right)\text{d}\alpha \text{d}r-cov\left(r_k,{\alpha }_i\right)-r_k^0{\alpha }_i^0}\right]}\\ -\frac{1}{2}{\displaystyle \underset{i=1}{\overset{TP}{\sum }}{\displaystyle \underset{j=1}{\overset{TP}{\sum }}{\lambda }_{ij}^{\left(22\right)}}\left[{\displaystyle \underset{D}{\int }{\alpha }_i{\alpha }_{j}p\left(\alpha ,r\right)\text{d}\alpha \text{d}r-\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right)-{\alpha }_i^0{\alpha }_j^0}\right]}.\end{array}$ (103)

In Equation (103), the quantities ${\lambda }_k^{\left(1\right)}$ , ${\lambda }_k^{\left(2\right)}$ , ${\lambda }_{k\mathcal{l}}^{\left(11\right)}$ , ${\lambda }_{k\mathcal{l}}^{\left(12\right)}$ , and ${\lambda }_{k\mathcal{l}}^{\left(22\right)}$ denote the respective Lagrange multipliers, and the factors 1/2 have been introduced for subsequent computational convenience.

Solving the equation $\partial H\left(p\right)/\partial p=0$ yields the following expression for the resulting MaxEnt distribution $p\left(z\right)$ :

$p_c\left(z\right)=\frac{1}{Z\left(b,\Lambda \right)}\exp \left(-b^{†}z-\frac{1}{2}{z}^{†}\Lambda z\right)$ , (104)

where the various vectors and matrices are defined as follows:

$z=\left(\begin{array}{c}r\\ \alpha \end{array}\right);b\triangleq \left(\begin{array}{c}{\lambda }^{\left(1\right)}\\ {\lambda }^{\left(2\right)}\end{array}\right);{\lambda }^{\left(1\right)}\triangleq \left(\begin{array}{c}{\lambda }_1^{\left(1\right)}\\ ·\\ {\lambda }_{TR}^{\left(1\right)}\end{array}\right);{\lambda }^{\left(2\right)}\triangleq \left(\begin{array}{c}{\lambda }_1^{\left(2\right)}\\ ·\\ {\lambda }_{TP}^{\left(1\right)}\end{array}\right);$ (105)

$\begin{array}{l}\Lambda \triangleq \left(\begin{array}{cc}{\Lambda }^{\left(11\right)}& {\Lambda }^{\left(12\right)}\\ {\left[{\Lambda }^{\left(12\right)}\right]}^{†}& {\Lambda }^{\left(22\right)}\end{array}\right);{\Lambda }^{\left(11\right)}\triangleq \left(\begin{array}{ccc}{\lambda }_{11}^{\left(11\right)}& ·& {\lambda }_{1,TR}^{\left(11\right)}\\ ·& {\lambda }_{k\mathcal{l}}^{\left(11\right)}& ·\\ {\lambda }_{TR,1}^{\left(11\right)}& ·& {\lambda }_{TR,TR}^{\left(11\right)}\end{array}\right);\\ {\Lambda }^{\left(12\right)}\triangleq \left(\begin{array}{ccc}{\lambda }_{11}^{\left(12\right)}& ·& {\lambda }_{1,TP}^{\left(12\right)}\\ ·& {\lambda }_{k\mathcal{l}}^{\left(12\right)}& ·\\ {\lambda }_{TR,1}^{\left(12\right)}& ·& {\lambda }_{TR,TP}^{\left(12\right)}\end{array}\right);{\Lambda }^{\left(22\right)}\triangleq \left(\begin{array}{ccc}{\lambda }_{11}^{\left(22\right)}& ·& {\lambda }_{1,TP}^{\left(22\right)}\\ ·& {\lambda }_{k\mathcal{l}}^{\left(22\right)}& ·\\ {\lambda }_{TP,1}^{\left(22\right)}& ·& {\lambda }_{TP,TP}^{\left(22\right)}\end{array}\right).\end{array}$ (106)

The normalization constant $Z\left(b,\Lambda \right)$ in Equation (104) is defined as follows:

$Z\left(b,\Lambda \right)={\displaystyle \underset{D}{\int }\exp \left(-b^{†}z-\frac{1}{2}{z}^{†}\Lambda z\right)\text{d}}z;\text{d}z\equiv \text{d}\alpha \text{d}r.$ (107)

In statistical mechanics, the normalization constant Z is called the partition function (or sum over states) and carries all of the information available about the possible states of the system, while the MaxEnt distribution $p\left(z\right)$ is called the canonical Boltzmann-Gibbs distribution. The integral in Equation (107) can be evaluated explicitly by conservatively extending the computational domain D to the entire multidimensional real vector space ${ℝ}^N$ , where $N\triangleq TR+TP$ , to obtain the following expression:

$Z_c\left(b,{\Lambda }_c\right)={\displaystyle \underset{{ℝ}^N}{\int }\exp \left(-b^{†}z-\frac{1}{2}{z}^{†}\Lambda z\right)\text{d}}z=\frac{{\left(2\pi \right)}^{N/2}}{\sqrt{Det\left(\Lambda \right)}}{\text{e}}^{\left(\frac{1}{2}{b}^{†}{\Lambda }^{-1}b\right)}.$ (108)

The Lagrange multipliers are determined in terms of the known information (means and covariances of parameters and responses) by differentiating the “free energy” $F\left(b,C\right)\triangleq -\ln Z\left(b,C\right)$ with respect to the components of the vector $b\triangleq {\left({\lambda }^{\left(1\right)},{\lambda }^{\left(2\right)}\right)}^{†}$ to obtain the following expressions:

$\frac{\partial F\left(b,C\right)}{\partial {\lambda }_k^{\left(1\right)}}=\frac{1}{Z}{\displaystyle \underset{D}{\int }{r}_k{\text{e}}^{\left(-b^{†}z-\frac{1}{2}{z}^{†}\Lambda z\right)}\text{d}z}={\displaystyle \underset{D_c}{\int }{r}_{k}p\left(\alpha ,r\right)\text{d}z}=r_k^0;k=1,\cdots ,TR;$ (109)

$\frac{\partial F\left(b,C\right)}{\partial {\lambda }_i^{\left(2\right)}}=\frac{1}{Z}{\displaystyle \underset{D}{\int }{\alpha }_i{\text{e}}^{\left(-b^{†}z-\frac{1}{2}{z}^{†}\Lambda z\right)}\text{d}z}={\displaystyle \underset{D}{\int }{\alpha }_{i}p\left(\alpha ,r\right)\text{d}\alpha \text{d}r}={\alpha }_i^0;i=1,\cdots ,TP.$ (110)

The results obtained in Equations (109) and (110) can be collectively written in vector-matrix form as follows:

$\frac{\partial F\left(b,C\right)}{\partial b}=z^0;z^0\triangleq {\left(r^0,{\alpha }^0\right)}^{†};r^0\triangleq \left(r_1^0,\cdots ,r_{TR}^0\right).$ (111)

On the other hand, it follows from Equation (108) that:

$\frac{\partial F\left(b,C\right)}{\partial b}=-\frac{\partial \left[\ln Z\left(b,C\right)\right]}{\partial b}=-{\Lambda }^{-1}b.$ (112)

The relations obtained in Equations (111) and (112) imply the following relation:

$b=-\Lambda z^0.$ (113)

Differentiating a second time the relation provided in Equation (109) or (110) yields the following relations:

$\begin{array}{c}\frac{{\partial }^{2}F\left(b,C\right)}{\partial {\lambda }_j^{\left(1\right)}\partial {\lambda }_k^{\left(1\right)}}=-\frac{1}{Z^2}\frac{\partial Z}{\partial {\lambda }_j^{\left(1\right)}}{r}_k^0+\frac{1}{Z}{\displaystyle \underset{D}{\int }\left(-r_j\right)r_k{\text{e}}^{\left(-b^{†}z-\frac{1}{2}{z}^{†}\Lambda z\right)}\text{d}z}\\ =-\mathrm{cov}\left(r_j,r_k\right);j,k=1,\cdots ,TR;\end{array}$ (114)

$\frac{{\partial }^{2}F\left(b,C\right)}{\partial {\lambda }_i^{\left(2\right)}\partial {\lambda }_k^{\left(1\right)}}=-\mathrm{cov}\left({\alpha }_i,r_k\right);i=1,\cdots ,TP;\text{}\text{}\text{}k=1,\cdots ,TR;$ (115)

$\frac{{\partial }^{2}F\left(b,C\right)}{\partial {\lambda }_i^{\left(2\right)}\partial {\lambda }_j^{\left(2\right)}}=-\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right);i,j=1,\cdots ,TP.$ (116)

The results obtained in Equations (114)-(116) can be collectively written in vector-matrix form as follows:

$\begin{array}{l}\frac{{\partial }^{2}F\left(b,C_c\right)}{\partial b\partial b}=-C;C\triangleq \left(\begin{array}{cc}{C}^{rr}& C^{r\alpha }\\ C^{\alpha r}& C^{\alpha \alpha }\end{array}\right);C^{rr}\triangleq {\left[\mathrm{cov}\left(r_j,r_k\right)\right]}_{TR\times TR};\\ C^{r\alpha }\triangleq {\left[\mathrm{cov}\left(r_k,{\alpha }_i\right)\right]}_{TR\times TP}={\left(C^{\alpha r}\right)}^{†};C^{\alpha \alpha }\triangleq {\left[\mathrm{cov}\left({\alpha }_i,{\alpha }_j\right)\right]}_{TP\times TP}.\end{array}$ (117)

On the other hand, it follows from Equation (112) that:

$\frac{{\partial }^{2}F\left(b,C\right)}{\partial b\partial b}=-{\Lambda }^{-1}.$ (118)

The relations obtained in Equations (117) and (118) imply the following relation:

${\Lambda }^{-1}=C.$ (119)

Introducing the results obtained in Equations (113) and (119) into Equation (104) and (108) yields the following expression for the MaxEnt distribution $p\left(z|z,C\right)$ :

$p\left(z|z,C\right)=\frac{{\left(2\pi \right)}^{-N/2}}{\sqrt{Det\left(C\right)}}\exp \left[-\frac{1}{2}{\left(z-z^0\right)}^{†}{C}^{-1}\left(z-z^0\right)\right].$ (120)

Appendix B: Construction of the “Constrained Minimum Discrepancy Distribution with a Given Function”

The aim is to obtain a mathematical expression for a distribution $p\left(z\right)$ , defined over a domain $z\in D$ , which is normalized to unity and satisfies K independent constraints of the following form:

${\displaystyle \underset{D}{\int }{g}_k\left(z\right)p\left(z\right)\text{d}}z=c_k<\infty ,k=1,\cdots ,K;\text{}\text{}\text{}$ (121)

with

${\displaystyle \underset{D}{\int }p\left(z\right)\text{d}}z=1$ . (122)

In addition, it is required that the normalized distribution $p\left(z\right)$ be “minimally discrepant from” (i.e., “fit as well as possible”) a known function $f\left(z\right)$ , which is also defined over the same domain $z\in D$ .

The discrepancy (i.e., “lack of fit”) between an approximate distribution $f\left(z\right)$ and a true distribution $p\left(z\right)$ , both defined over the same domain $z\in D$ , is denoted below as $\delta \left(f|p\right)$ , and is defined [see, e.g., ] as follows:

$\delta \left(f|p\right)\triangleq {\displaystyle \underset{D}{\int }p\left(z\right)\log \frac{p\left(z\right)}{f\left(z\right)}\text{d}}z.$ (123)

The discrepancy $\delta \left(f|p\right)$ is also called, see e.g., [22] , the “Kulback-Leibler divergence or relative entropy.” The normalized distribution $p\left(z\right)$ which minimizes the discrepancy defined in Equation (123) while satisfying the constraints expressed by Equations (121) and (122) is obtained by minimizing the Lagrangian functional $F\left(p\right)$ defined below:

$F\left(p\right)\triangleq {\displaystyle \underset{D}{\int }p\left(z\right)\log \frac{p\left(z\right)}{f\left(z\right)}\text{d}}z+{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k}\left[{\displaystyle \underset{D}{\int }{g}_k\left(z\right)p\left(z\right)\text{d}}z-c_k\right]+c\left[{\displaystyle \underset{D}{\int }p\left(z\right)\text{d}}z-1\right].$ (124)

In Equation (124), the quantities c and ${\theta }_1,\cdots ,{\theta }_k,\cdots ,{\theta }_{TR}$ are Lagrange multipliers which remain to be determined. The stationary value of $F\left(p\right)$ is provided by the probability $p\left(z\right)$ where the first variation of $F\left(p\right)$ vanishes, i.e., $p\left(z\right)$ is the solution of the following equation

$\begin{array}{c}\delta F\left(p\right)\triangleq \frac{\text{d}}{\text{d}\epsilon }{\left\{F\left[p\left(z\right)+\epsilon \tau \left(z\right)\right]\right\}}_{\epsilon =0}\\ ={\displaystyle \underset{D}{\int }\left[\log \frac{p\left(z\right)}{f\left(z\right)}+{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k{g}_k\left(z\right)+\left(c+1\right)}\right]\tau \left(z\right)\text{d}}z=0\text{}\text{}\text{}.\text{}\text{}\end{array}$ (125)

The desired posterior distribution is the normalized solution of Equation (125), having the following expression:

$p\left(z\right)=\frac{1}{Z}f\left(z\right)\exp \left[-{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k{g}_k\left(z\right)}\right];Z\triangleq {\displaystyle \underset{D}{\int }f\left(z\right)\exp \left[-{\displaystyle \underset{k=1}{\overset{TR}{\sum }}{\theta }_k{g}_k\left(z\right)}\right]\text{d}}z\text{}.\text{}$ (126)

It is important to note that in the limiting case when the function $f\left(z\right)$ is the uniform distribution $f\left(z\right)=\text{constant}$ , the distribution $p\left(z\right)$ becomes a MaxEnt-distribution. In other words, the distribution which is least discrepant from (i.e., “best-fits”) the uniform distribution is the MaxEnt distribution.

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

References

[1]	Sasaki, Y.K. (1955) A Fundamental Study of the Numerical Prediction Based on the Variational Principle. Journal of the Meteorological Society of Japan, 33, 262-275. https://doi.org/10.2151/jmsj1923.33.6_262
[2]	Sasaki, Y.K. (1958) An Objective Analysis Based on the Variational Method. Journal of the Meteorological Society of Japan, 36, 77-88. https://doi.org/10.2151/jmsj1923.36.3_77
[3]	Eckart, C. (1960) Hydrodynamics of Ocean and Atmosphere. Pergamon Press, Oxford.
[4]	Humi, I., Wagschal, J.J. and Yeivin, Y. (1964) Multi-Group Constants from Integral Data. Proceedings of the 3rd International Conference on the Peaceful Uses of Atomic Energy, Geneva, 31 August-9 September 1964, 398.
[5]	Lewis, J.M., Lakshmivarahan, S. and Dhall, S.K. (2006) Dynamic Data Assimilation: A Least Square Approach. Cambridge University Press, Cambridge. https://doi.org/10.1017/CBO9780511526480
[6]	Lahoz, W., Khattatov, B. and Ménard, R. (2010) Data Assimilation: Making Sense of Observations. Springer Verlag, Heidelberg. https://doi.org/10.1007/978-3-540-74703-1
[7]	Faragó, I., Havasi, á. and Zlatev, Z. (2013) Advanced Numerical Methods for Complex Environmental Models: Needs and Availability. Bentham Science Publishers, Bussum. https://doi.org/10.2174/97816080577881130101
[8]	Cacuci, D.G., Navon, M.I. and Ionescu-Bujor, M. (2014) Computational Methods for Data Evaluation and Assimilation. Chapman & Hall/CRC, Boca Raton.
[9]	Rowlands, J., et al. (1973) The Production and Performance of the Adjusted Cross-Section Set FGL5. Proceedings International Symposium on Physics of Fast Reactors, Tokio. http://www.cern.ch/root
[10]	Gandini, A. and Petilli, M. (1973) A Code Using the Lagrange’s Multipliers Method for Nuclear Data Adjustment. CNEN-RT/FI(73)39, Comitato Nazionale Energia Nucleare, Rome.
[11]	Kuroi, H. and Mitani, H. (1975) Adjustment to Cross-Section Data to Fit Integral Experiments by Least Squares Method. Journal of Nuclear Science and Technology, 12, 663-680. https://doi.org/10.1080/18811248.1975.9733172
[12]	Dragt, J.B., Dekker, J.W.M., Gruppelaar, H. and Janssen, A.J. (1977) Methods of Adjustment and Error Evaluation of Neutron Capture Cross Sections. Nuclear Science and Engineering, 62, 117-129. https://doi.org/10.13182/NSE77-3
[13]	Cacuci, D.G. (2014) Predictive Modeling of Coupled Multi-Physics Systems: I. Theory. Annals of Nuclear Energy, 70, 266-278. https://doi.org/10.1016/j.anucene.2013.11.027
[14]	Cacuci, D.G. (2019) BERRU Predictive Modeling: Best-Estimate Results with Reduced Uncertainties. Springer, Berlin. https://doi.org/10.1007/978-3-662-58395-1
[15]	Jaynes, E.T. (1957) Information Theory and Statistical Mechanics. Physical Review, 106, 620-630. https://doi.org/10.1103/PhysRev.106.620
[16]	Cacuci, D.G. (1981) Sensitivity Theory for Nonlinear Systems: I. Nonlinear Functional Analysis Approach. Journal of Mathematical Physics, 22, 2794-2812. https://doi.org/10.1063/1.525186
[17]	Cacuci, D.G. (2016) The Second-Order Adjoint Sensitivity Analysis Methodology for Nonlinear Systems—I: Theory. Nuclear Science and Engineering, 184, 16-30. https://doi.org/10.13182/NSE16-16
[18]	Cacuci, D.G. (2023) Perspective on Predictive Modeling: Current Status, New High-Order Methodology and Outlook for Energy Systems. Energies, 16, Article 933. https://doi.org/10.3390/en16020933
[19]	Cacuci, D.G. (2023) Second-Order MaxEnt Predictive Modeling Methodology. II: Probabilistically Incorporated Computational Model Responses (2nd-BERRU-PMP). American Journal of Computational Mathematics, 13, No. 2.
[20]	Shun, Z. and McCullagh, P. (1995) Laplace Approximation of High Dimensional Integrals. Journal of the Royal Statistical Society: Series B (Methodological), 57, 749-760. https://doi.org/10.1111/j.2517-6161.1995.tb02060.x
[21]	Shannon, C.E. (1948) A Mathematical Theory of Communication. The Bell System Technical Journal, 27, 379-423. https://doi.org/10.1002/j.1538-7305.1948.tb01338.x
[22]	Bernardo, J.M. and Smith, A.F.M. (2006) Bayesian Theory. John Wiley & Sons, Chichester.