Uncertainty Quantification in Case of Imperfect Models: A Review  *** Running
title: Uncertainty quantification: a review

Sebastian Kersting††† Corresponding author. Tel: +49-6151-16-23374, Fax:
+49-6151-16-23381 and Michael Kohler

Fachbereich Mathematik, Technische Universität Darmstadt, Schlossgartenstr. 7,
64289 Darmstadt, Germany, email<EMAIL_ADDRESS><EMAIL_ADDRESS>
December 11, 2020

Abstract
Uncertainty quantification of complex technical systems is often based on a
computer model of the system. As all models such a computer model is always
wrong in the sense that it does not describe the reality perfectly. The
purpose of this article is to give a review of techniques which use observed
values of the technical systems in order to take into account the inadequacy
of a computer model in uncertainty quantification. The techniques reviewed in
this article are illustrated and compared by applying them to applications in
mechanical engineering.

Key words and phrases: Estimated input distributions, imperfect models,
improved surrogate models, uncertainty quantification.

##  1 Introduction

Uncertainty quantification is a major research field of statistical methods
with applications in engineering sciences. Methods of uncertainty
quantification are often used to analyse experiments with technical systems.
These experiments can be described by $\mathbb{R}^{d}\times\mathbb{R}$-valued
(random) variables $(X,Y)$, where $X$ describes input parameters of the
experiments and $Y$ describes the outcome of the experiment. E.g. if one
considers drop tests of a spring damper as described in Subsection 2.1 the
experimental outcome depends on the drop height, i.e. in this case the input
dimension $d$ is equal to $1$ and the measured outcome is the maximal relative
compression.

Usually to conduct real experiments is expensive and time consuming. To
circumvent this problem computer models that simulate the experiment with the
technical system are playing crucial role. Santner et al., (2018) and Fang et
al., (2006) provide an overview of methods for the design and analysis of
experiments conducted with these computer models.

Computer models are often computationally expensive and thus it is not
possible to conduct a large quantity of computer experiments. To circumvent
this problem one can use so-called surrogate models. There exists different
methods to estimate a surrogate model, e.g. Bucher and Bourgund, (1990), Kim
and Na, (1997). Das and Zheng, (2000) used quadratic response surfaces.
Hurtado, (2004), Deheeger and Lemaire, (2010), Bourinet et al., (2011)
investigated surrogate models in context of support vector machines,
Papadrakakis and Lagaros, (2002) concentrated on neural networks, Kaymaz,
(2005) and Bichon et al., (2008) used kriging. Usually these surrogate models
do not take an inherent error of the computer model into account and thus are
themself imperfect in case of an imperfect computer model.

There is a number of different objectives in methods of uncertainty
quantification. The main objectives are computer model calibration, computer
model validation and quantifying the results of an experiment. Kennedy and
O’Hagan, (2001) and Han et al., (2009) proposed a Bayesian method for computer
model calibration and model validation based on Gaussian processes, Higdon et
al., (2013) proposed a similar approach as Kennedy and O’Hagan, (2001), but
used the ensemble Kalman filter. Plumlee, (2017) also introduced a Bayesian
computer model calibration method similar to the approach of Kennedy and
O’Hagan, (2001), but their prior distribution on the bias (model discrepancy)
is modeled orthogonal to the gradient of the computer model. Bayarri et al.,
(2007) proposed a validation method in case of time dependent systems, where
they model the discrepancy between the computer experiments and the outcome of
the technical system by a Gaussian process. Wang et al., (2009) proposed a
Bayesian method for model validation. Gu and Wang, (2018) proposed the so-
called scaled Gaussian Process and use it for model calibration. They claim,
that the method bridges the gap between the $L_{2}$ calibration and the
Gaussian process calibration. Damblin et al., (2018) proposed a method for
computer model calibration based on perfect and imperfect computer models. Tuo
and Wu, (2015) proposed the $L_{2}$ calibration model, to identify unknown
computer model parameter. They also pointed out that the Bayesian approach of
Kennedy and O’Hagan, (2001) might fail in case of an imperfect computer model,
since there exists no values of the parameters which fit the technical system
perfectly. Plumlee, (2019) proposed a frequentist method for computer model
calibration by constructing confidence intervals for the model parameter.

Based on perfect computer models, Devroye et al., (2013) and Bott et al.,
(2015) derived consistency and rate of convergence results for density
estimators based on surrogate models, Felber et al., 2015a  and Felber et al.,
2015b  proposed a method for the adaptive choice of smoothing parameters for
density estimators based on surrogate models and Kohler and Krzyżak, (2019)
proposed a method to estimate quantiles based on surrogate models.

The articles mentioned above either assume that the underlying computer model
is perfect, or only deal with computer model calibration. Normally, computer
models are imperfect, i.e. they do not predict the outcome of real experiments
perfectly, e.g. because of a relaxation of underlying physical dependencies
(typically neglecting the friction or by considering it to be constant) to
reduce complexity or because of missing knowledge about the technical system.
In this article we give an overview of existing methods based on imperfect
computer models with the following goals. The first aim is to construct an
improved surrogate model and based on that quantify the outcome of the
experiment $Y$ either by density or quantile estimation. The second is to
quantify the model error and thus enable the comparison of different computer
models. The last aim is to quantify the influence of the computer model error
on the quantification of $Y$.

###  1.1 Notation

Throughout this paper we use the following notation: $\mathbb{N}$,
$\mathbb{R}$ and $\mathbb{R}_{+}$ are the sets of positive integers, real
numbers, and nonnegative real numbers, respectively. For $z\in\mathbb{R}$ we
denote the smallest integer greater than or equal to $z$ by $\lceil z\rceil$.
For $x\in\mathbb{R}^{d}$ we denote the $i$-th component of $x$ by $x^{(i)}$.
If $X$ is a random variable, then ${\mathbf{P}}_{X}$ is the corresponding
distribution, i.e., the measure associated with the random variable. Let
$D\subseteq\mathbb{R}^{d}$ and let $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ be
a real-valued function defined on $\mathbb{R}^{d}$. We write $x=\arg\min_{z\in
D}f(z)$ if $\min_{z\in{\mathcal{D}}}f(z)$ exists and if $x$ satisfies

$x\in D\quad\mbox{and}\quad f(x)=\min_{z\in{\mathcal{D}}}f(z).$

If $A$ and $B$ are sets, with $A\subseteq B$. then $I_{A}\colon B\to\\{0,1\\}$
is the indicator function corresponding to $A$, i.e. the function which takes
on the value 1 on A and is 0 elsewhere.

###  1.2 Outline

The outline of this paper is as follows: In Section 2 we introduce two
technical systems, which will be used in the following sections to illustrate
and discuss the described methods. In Section 3 we describe typical data
models and a method used to generate additional input values. In Section 4
uncertainty quantification methods based on improved surrogate models are
described and illustrated. In Section 5 methods used to compare different
computer simulations are described and illustrated. In Section 6 two methods
which can be used to quantify the influence of the computer model error are
described and illustrated.

##  2 Technical systems

In the following, two technical systems are introduced. The data generated by
experiments with these systems will be used to illustrate the methods of
uncertainty quantification.

###  2.1 Modular active spring damper system

The modular active spring damper system (German acronym MAFDS) shown in Figure
1 is a suspension system that was designed with similar specifications and
requirements as an air plane landing gear, although it is not a landing gear
substitute. It was developed in the collaborative research centre SFB 805 at
the Technische Universität Darmstadt in order to investigate uncertainty in a
load-bearing structural system when predicting the dynamic response.

Figure 1:  A photo of the MAFDS and its experimental test setup.

Its main components are an upper truss structure, a lower truss structure,
guidance links that enable relative translation of the truss structures in
vertical direction and a spring-damper. A detailed system description can be
found in Feldmann and Platz, (2019) and Mallapur and Platz, (2019). Dynamic
drop tests can be carried out similar to landing gear testing: the MAFDS is
lifted up and dropped by a variable drop height $h$. In a total of $100$
experiments, the inputs $h$ were chosen normally distributed with mean
$0.05~{}\mathrm{m}$ and a standard deviation of $0.0057~{}\mathrm{m}$. As
system output we regard the maximal relative compression $z_{\mathrm{r,max}}$
between the upper and the lower truss during the drop test.

Modelling of the system at hand yielded two different computer models to
compute the maximum relative compression $z_{\mathrm{r,max}}$ to a
corresponding input value. With both computer models we conduct $500$
simulations each and furthermore generate $16.000.000$ additional input values
that will be subsequently used in this paper.

###  2.2 Piezo-elastic beam

The second technical system under consideration is a lateral vibration
attenuation system with piezo-elastic supports shown in Figure 2.

Figure 2:  A CAD model of the lateral vibration attenuation system with piezo-
elastic supports and a sectional view of one of the piezo-elastic supports,
cf. Li et al., (2017).

This system consists of a beam with circular cross-section embedded in two
piezo-elastic supports A and B where support A is used for lateral beam
vibration excitation and support B is used for lateral beam vibration
attenuation, as proposed in Götz et al., (2016) and Götz, (2019). Both
supports are located at opposite beam ends and each consists of one elastic
membrane-like spring element made of spring steel, two piezoelectric stack
transducers arranged orthogonally to each other and mechanically pre-stressed
with disc springs as well as the relatively stiff axial extension made of
hardened steel that connects the piezoelectric transducers with the beam. For
vibration attenuation in support B, optimally tuned electrical shunt circuits
are connected to the piezoelectric transducers.

As system output serves the maximum lateral vibration amplitude $y$ in an
experiment with this system. Five quantities of the membrane in the
attenuation system vary during the manufacturing of the piezo-elastic supports
and influence the maximal vibration amplitude: the lateral stiffness in
$y$-direction $k_{lat,y}$ and in $z$-direction $k_{lat,z}$, the rotatory
stiffness in $y$-direction $k_{rot,y}$ and in $z$-direction $k_{rot,z}$, and
the height of the membrane $h_{x}$. Ten measurements of corresponding
parameters together with the experimental outcome $y$ for ten built systems
are available, which are listed in Table 1.

| 1 | 2 | 3 | 4 | 5 | 6 | 7 | 8 | 9 | 10
---|---|---|---|---|---|---|---|---|---|---
$k_{rot,y}\times 10^{2}$ | 1.31 | 1.34 | 1.31 | 1.23 | 1.14 | 1.29 | 1.35 | 1.28 | 1.04 | 1.20
$k_{rot,z}\times 10^{2}$ | 1.31 | 1.28 | 1.43 | 1.25 | 1.30 | 1.34 | 1.22 | 1.16 | 1.18 | 1.11
$k_{lat,y}\times 10^{7}$ | 3.27 | 3.28 | 3.35 | 3.29 | 3.22 | 3.26 | 3.19 | 3.54 | 3.21 | 3.42
$k_{lat,z}\times 10^{7}$ | 3.07 | 3.22 | 3.29 | 3.25 | 3.30 | 3.18 | 3.16 | 3.51 | 3.37 | 3.44
$h_{x}\times 10^{-4}$ | 6.79 | 6.77 | 6.82 | 6.80 | 6.79 | 6.76 | 6.81 | 6.74 | 6.68 | 6.84
$y\times 10^{1}$ | 1.45 | 1.42 | 1.44 | 1.42 | 1.43 | 1.35 | 1.47 | 1.32 | 1.31 | 1.63
Table 1: Measurement data for ten built systems. The values of $k_{rot,y}$ and
$k_{rot,z}$ are given in $[Nm/\operatorname{rad}]$, the values of $k_{lat,y}$
and $k_{lat,z}$ are given in $[N/m]$, the values of $h_{x}$ are given in $[m]$
and the values of the maximal vibration amplitude $y$ are given in
$[\frac{m}{s^{2}}/V]$.

A computer model is available with which we can compute the maximal vibration
amplitude $y$ to the corresponding input values of $k_{rot,y}$, $k_{rot,z}$,
$k_{lat,y}$, $k_{lat,z}$ and $h_{x}$. Apart from the input values displayed in
Table 1 additional input values are generated. Therefore, we assume that the
input values are independent multivariate normally distributed and we generate
a sample of additional input parameters as described in Subsection 3.1. For
simulations of the computer model, a total of $200$ additional input values
are generated. In addition, $16.000.000$ additional input values are generated
that will be subsequently used in this paper.

##  3 Data models

In an application, usually one of the following two different data models are
available:

  1. (i)

Full information model: In this case, a computer simulation
$m\colon\mathbb{R}^{d}\to\mathbb{R}$ of the experiment with the technical
system is available and can be evaluated at arbitrary values, an experimental
data set

$(X_{1},Y_{1}),\ldots,(X_{n},Y_{n})$ (1)

of size $n\in\mathbb{N}$ is available and the underlying distribution
${\mathbf{P}}_{X}$ of $X$ is known thus one can generate an additional data
set

$X_{n+1},\ldots,X_{n+L_{n}+N_{1,n}+N_{1,n}}$ (2)

of additional $L_{n}+N_{1,n}+N_{2,n}$ input values.

  2. (ii)

Experimental setup model: In this case, a computer simulation
$m\colon\mathbb{R}^{d}\to\mathbb{R}$ of the experiment with the technical
system and the experimental data set (1) of size $n\in\mathbb{N}$ are
available. Here the underlying distribution ${\mathbf{P}}_{X}$ of $X$ is
unknown.

In the following section we will describe a method to generate an additional
sample of input values. Thus, to simplify the presentation we will assume that
the full information model is available for the description of the methods.
Note that for the MAFDS, cf. Subsection 2.1, data model (i) is available and
for the piezo-beam, cf. Subsection 2.2, data model (ii) is available.

###  3.1 Generating additional input values

To apply methods of uncertainty quantification, often a large sample of input
data is necessary. In some settings (especially if the input data is random as
well as the outputs) only a small data set of input values is available. In
this case it is necessary to generate more input data. But this is sometimes
nearly as expensive and time consuming as producing more output data, e.g. in
the application from Subsection 2.2 one has to produce a new disk spring and
measure the four stiffness parameters and the height parameter. Thus
generating more experimental input data is often not possible.

The randomness in input data is often induced by production margins or
measurement errors and thus making a realistic assumption on the underlying
class of distributions is often possible. A variety of methods for data
generation in case of known distributions is given in Devroye, (1986). In the
following, the in Kersting and Kohler, (2019) proposed method to generate
multivariate normally distributed random variables is described:

Assume $X$ is a multivariate normally distributed random variable with values
in $\mathbb{R}^{d}$. Furthermore assume that an independent and identically
distributed sample

$X_{1},\ldots,X_{n}$ (3)

of $X$ is available.

Then estimate the parameters of the distribution of $X$ by a maximum
likelihood estimate defined by

$\hat{\mu}=\frac{1}{n}\sum_{i=1}^{n}X_{i}$ (4)

and

$\hat{\Sigma}=\left(\frac{1}{n}\sum_{k=1}^{n}(X_{k}^{(i)}-\hat{\mu}^{(i)})(X_{k}^{(j)}-\hat{\mu}^{(j)})\right)_{1\leq
i,j\leq d}.$ (5)

The estimated values are then treated as if they were the real values. Then
one can generate a sample of size $N_{n}\in\mathbb{N}$ which is independent
and multivariate normally distributed with mean $\hat{\mu}$ and covariance
$\hat{\Sigma}$ by firstly generating an independent sample
$Z_{1},\ldots,Z_{N_{n}}$ of $d$-dimensional vectors, where for each vector the
components are independent and standard normally distributed, and set for
every $i=1,\ldots,N_{n}$

$\bar{X}_{i}=\hat{O}\hat{\Lambda}^{1/2}Z_{i}+\hat{\mu},$ (6)

where $\hat{O}$ and $\hat{\Lambda}$ are defined by the eigendecomposition

$\hat{\Sigma}=\hat{O}\hat{\Lambda}\hat{O}^{T}$

of $\hat{\Sigma}$. Here
$\hat{\Lambda}=\operatorname{diag}(\hat{\lambda}_{1},\ldots,\hat{\lambda}_{d})$
is a diagonal matrix consisting of eigenvalues of $\hat{\Sigma}$ and $\hat{O}$
is a orthogonal matrix whose columns are eigenvectors of $\hat{\Sigma}$.

##  4 Uncertainty quantification based on improved surrogate models

As discussed above, computer models are usually imperfect. In this section, a
method to estimate an improved surrogate model and it’s application in
uncertainty quantification is described.

###  4.1 Estimating an improved surrogate model

Computer simulations of experiments with a technical system are often complex
and computationally expensive to evaluate. But for a proper analysis of the
underlying technical system it is necessary to generate a large sample of
these computer experiments. In addition, the inherent error in the computer
simulation, cf. Section 1, will bias the results. To circumvent these
problems, a solution is to compute an improved surrogate model and use it
instead of the computer simulation.

The first step is to generate a sample of input values evaluated with the
computer simulation. Kohler and Krzyżak, (2017) and Götz et al., (2018) used a
sample of the input quantity $X$, Wong et al., (2017) suggested to use data
generated in a specific range by Latin Hypercube sampling (see McKay et al.,
(1979)), if the input data is not random, and Kersting and Kohler, (2019) used
a sample uniformly distributed on a centered cube to construct an estimator of
a computer model.

The second step is to estimate a surrogate model of the computer simulation.
There exists a vast variety of methods, some already mentioned above. One
possibility is to use (penalized) least squares estimates, defined by

$\hat{m}_{L_{n}}(\cdot)=\arg\min_{f\in{\cal
F}}\frac{1}{L_{n}}\sum_{i=1}^{L_{n}}|f(X_{i})-m(X_{i})|^{2}+pen_{n}^{2}(f),$
(7)

where ${\cal F}$ is a set of functions,
$(X_{1},m(X_{1})),\ldots,(X_{L_{n}},m(X_{L_{n}}))$ is a set of input values
evaluated with the computer model of size $L_{n}\in\mathbb{N}$ and
$pen_{n}^{2}(\cdot)$ is a penalty term which usually penalizes the roughness
of the function and which is positive for each $f\in{\cal F}$, i.e.
$pen_{n}^{2}(f)\geq 0$. If the input dimension is smaller or equal to $3$ than
for example smoothing spline estimates can be applied as shown in Kohler and
Krzyżak, (2017). For bigger input dimensions neural network estimates can be
used as in Götz et al., (2018) and Kersting and Kohler, (2019).

As discussed, usually every computer model has an inherent error. To improve
the surrogate model an estimator of the residuals can be constructed by first
calculating the residuals of the surrogate model on the experimental data

$\epsilon_{i}=Y_{i}-\hat{m}_{L_{n}}(X_{i})\quad(i=1,\ldots,n)$ (8)

and by applying then a (penalized) least squares estimate to this sample, i.e.
by computing

$\hat{m}_{n}^{\epsilon}(\cdot)=\arg\min_{f\in\bar{{\cal
F}}}\frac{1}{n}\sum_{i=1}^{n}|f(X_{i})-\epsilon_{i}|^{2}+pen_{n}^{2}(f),$ (9)

where $\bar{{\cal F}}$ is a set of functions and $pen_{n}^{2}(\cdot)$ is a
penalty term for which $pen_{n}^{2}(f)\geq 0$ holds, for every $f\in\bar{{\cal
F}}$.

In the last step, the improved surrogate model is defined by

$\hat{m}_{n}(x)=\hat{m}_{L_{n}}(x)+\hat{m}_{n}^{\epsilon}(x)\quad(x\in\mathbb{R}).$
(10)

In case that only a small sample of experimental data is available, i.e. $n$
is small, the residual estimator usually does not yield sufficient results.
Then one can use a weighted (penalized) least squares estimate instead,
defined by

$\hat{m}_{n}^{\epsilon}(\cdot)=\arg\min_{f\in\bar{{\cal
F}}}\left(\frac{w^{(n)}}{n}\sum_{i=1}^{n}|f(X_{i})-\epsilon_{i}|^{2}+\frac{(1-w^{(n)})}{N_{1,n}}\sum_{i=1}^{N_{1,n}}|f(X_{i})-0|^{2}+pen_{n}^{2}(f)\right),$
(11)

where $w^{(n)}\in[0,1]$ is a weighting term, which should be chosen data
dependent and $X_{n+L_{n}+1},\ldots,X_{n+L_{n}+N_{1,n}}$ is a set of
additional input values of size $N_{1,n}\in\mathbb{N}$. Note that in (7) and
(9) different sets of functions are used. In an application this is usually
due to different chosen function class parameters.

###  4.2 Density estimation based on an improved surrogate model

One approach to quantifying the uncertainty in the outcome of an experiment
with a technical system is to estimate the underlying density of $Y$. If the
density $g\colon\mathbb{R}\to\mathbb{R}$ is known, we can calculate the
probability ${\mathbf{P}}\\{Y\in B\\}$ for any given set
$B\subseteq\mathbb{R}$ by

$\int_{B}g(x)\;dx.$ (12)

In the following section we will describe a method to construct an estimator
$\hat{g}$ of the density $g$ of $Y$. A simple approach is to apply the kernel
density estimate of Rosenblatt, (1956) and Parzen, (1962) defined by

$g(y)=\frac{1}{n\cdot
h_{n}}\cdot\sum_{i=1}^{n}K\left(\frac{y-Y_{i}}{h_{n}}\right).$ (13)

Usually in an application the sample size of experimental data will be too
small to achieve satisfying results. Instead one can use an improved surrogate
model as described in Subsection 4.1.

Kohler and Krzyżak, (2017), Götz et al., (2018) and Kersting and Kohler,
(2019) proposed similar methods to construct a density estimate based on an
improved surrogate model, where the main difference is that the first two
assumed that the full information model is available, where the last assumed
that the experimental setup model is given and that Kohler and Krzyżak, (2017)
used smoothing spline estimates, where Götz et al., (2018) and Kersting and
Kohler, (2019) used least squares neural network estimates to construct the
improved surrogate model. Summarizing, the method works as follows:

Estimate an improved surrogate model as discussed in Subsection 4.1. Next
apply a kernel density estimator to an independent and identically distributed
sample of $(X,\hat{m}_{n}(X))$ of size $N_{2,n}\in\mathbb{N}$, i.e. to

$\hat{g}_{N_{2,n}}(y)=\frac{1}{N_{2,n}\cdot
h_{N_{2,n}}}\cdot\sum_{i=1}^{N_{2,n}}K\left(\frac{y-\hat{m}_{n}(X_{n+L_{n}+N_{1,n}+i})}{h_{N_{2,n}}}\right),$
(14)

where $h_{N_{2,n}}>0$ so-called bandwidth and $K\colon\mathbb{R}\to\mathbb{R}$
the kernel function are parameters of the estimate. To choose the parameter
$h_{N_{2,n}}$, a data dependent approach can be used as shown in Felber et
al., 2015a  or Felber et al., 2015b .

###  4.3 Quantile estimation based on an improved surrogate model

Another approach to quantify the uncertainty in the outcome of an experiment
with a technical system is to estimate its $\alpha$-quantile, defined by

$q_{Y,\alpha}=\inf\left\\{y\in\mathbb{R}:G(y)\geq\alpha\right\\},$ (15)

where $\alpha\in(0,1)$ and $G$ is the cumulative distribution function of $Y$
given by

$G(y)={\mathbf{P}}\left\\{Y\leq y\right\\}.$ (16)

A simple and straight-forward approach is to use a Monte Carlo estimator. Enss
et al., (2016) and Kohler and Krzyżak, (2018) proposed to use an estimator
based on a surrogate model, where they assumed that the computer simulation
fits the reality perfectly. As discussed before, a computer simulation is
almost always imperfect and thus a simple modification is to use an improved
surrogate model. The estimator is then constructed as follows:

Estimate an improved surrogate model as discussed in Subsection 4.1. Next
apply a Monte Carlo quantile estimator to an independent and identically
distributed sample of $(X,\hat{m}_{n}(X))$ of size $N_{2,n}\in\mathbb{N}$,
defined by

$\hat{q}_{\hat{m}_{n}(X),N_{2,n},\alpha}=\inf\left\\{y\in\mathbb{R}:\hat{G}_{\hat{m}_{n}(X),N_{2,n}}(y)\geq\alpha\right\\},$
(17)

where

$\hat{G}_{\hat{m}_{n}(X),N_{2,n}}(y)=\frac{1}{N_{2,n}}\sum_{i=1}^{N_{2,n}}I_{(-\infty,m_{n}(X_{n+L_{n}+N_{1,n}+i})]}(y).$
(18)

Note that this is equivalent to choosing the $\lceil N_{2,n}\cdot\alpha\rceil$
biggest value of $\hat{m}_{n}(X_{n+L_{n}+N_{1,n}+i})$,
$\ldots$,$\hat{m}_{n}(X_{n+L_{n}+N_{1,n}+N_{2,n}})$.

###  4.4 Application

In the following, we will illustrate and discuss the above described methods
by applying them to data of the technical systems described in Section 2. In
case of the MAFDS we restrict ourselves to use only $10$ randomly chosen
experimental data points of the available $100$.

####  4.4.1 Estimation of an improved surrogate model

Many of the in the following described methods of uncertainty quantification
need an (improved) surrogate model of the computer model. If not stated
otherwise the (improved) surrogate model will be estimated as follows: In case
of the MAFDS, cf. Subsection 2.1, the input dimension $d$ is $1$. Here we use
in (7) a smoothing spline estimator as implemented by the MATLAB method csaps,
where the smoothing parameter is chosen by generalized cross-validation. In
(10) the smoothing parameter and the weighting parameter
$w^{(n)}\in\\{0,0.1,0.2,\ldots,1\\}$ are chosen simultaneously by a $5$-fold
cross validation, where the empirical $L_{2}$ risk is calculated only on the
experimental data. In case of the piezo-beam, cf. Subsection 2.2, we use fully
connected feed forward neural networks from the MATLAB Deep Learning Toolbox
to estimate the surrogate model. The network topology is chosen data dependent
by a splitting of the sample, where we use $2/3$ training and $1/3$ testing
data. We use networks with $\\{1,5,10\\}$ hidden layer which have
$\\{1,6,11,16,21\\}$ neurons per layer. To estimate an improved surrogate
model we use a special class of sparsely connected neural networks. A detailed
description of these neural networks and their implementation can be found in
Kersting and Kohler, (2019).

For the MAFDS the resulting improved surrogate models, together with the
experimental data and the surrogate models based on the two computer models
are displayed in the left plots of Figures 3 and 4. The method estimates
(nearly) identical surrogate models for both computer simulation. The
resulting surrogate model fits the experimental data better than both computer
simulations.

Figure 3: Left plot contains $10$ randomly chosen experimental data points of
the MAFDS (circles) used to construct the improved surrogate model, the
surrogate model based on computer model 1 (dashed line) and the improved
surrogate model (solid line). Right plot contains a kernel density estimator
based on the $10$ experimental data (gray solid line), a kernel density
estimator based on the computer model 1 (black dashed line) and kernel density
estimator based on an improved surrogate model (black solid line). The $10$
experimental data points are indicated on the x-axis. Furthermore, for
comparison a kernel density estimator based on all available experimental data
(gray dashed line).

Figure 4: Left plot contains $10$ randomly chosen experimental data points of
the MAFDS (circles) used to construct the improved surrogate model, the
surrogate model based on computer model 2 (dashed line) and the improved
surrogate model (solid line). Right plot contains a kernel density estimator
based on the $10$ experimental data (gray solid line), a kernel density
estimator based on the computer model 2 (black dashed line) and kernel density
estimator based on an improved surrogate model (black solid line). The $10$
experimental data points are indicated on the x-axis. Furthermore, for
comparison a kernel density estimator based on all available experimental data
(gray dashed line).

Since the piezo-beam depends on $5$ input parameters, we are not able to
visualize the result for the piezo-beam.

####  4.4.2 Density estimation based on an improved surrogate model

We apply the density estimation method to data of both technical systems. The
surrogate models are constructed as described in Subsection 4.1. The results
are displayed in Figures 3, 4 and 5.

Figure 5: Plot compares density estimator of piezo-beam, cf. Subsection 2.2. A
kernel density estimator based on all available experimental data (gray solid
line), a kernel density estimator based on the computer model (black dashed
line) and kernel density estimator based on an improved surrogate model (black
solid line). The values of the outcome of the $10$ experiments are indicated
on the x-axis.

In case of the MAFDS, the resulting improved surrogate models are (almost)
equal for both computer models, although the estimators based solely on the
computer models differ a lot. For both computer models, the improved surrogate
model seems to correct their prediction error and yields a more accurate
estimation of the real density. In case of the piezo-beam, the resulting
improved surrogate model predict higher values than the computer simulation,
which fits the experimental data better.

To conclude, in both applications the improved surrogate model is able to
correct the computer model error and yields a better prediction of the real
experiments, resulting in an density estimation which fits the reality better
than the density estimation based on computer simulations, or the density
estimator based solely on experimental data.

####  4.4.3 Quantile estimation based on an improved surrogate model

We apply the quantile estimation method described in Subsection 4.3 to the
data of both technical systems. The surrogate models are constructed as
described in Subsection 4.1. In case of the MAFDS we restrict ourselves to use
only $10$ of the available $100$ experimental data points to estimate the
improved surrogate model. For both technical systems we are interested in the
$0.95$-quantile. The resulting estimates are shown in Table 2.

| MAFDS | piezo-beam
---|---|---
computer model 1 | $0.0849$ | $14.6994$
computer model 2 | $0.0849$ |
Table 2:  Results of the Monte Carlo estimator based on an improved surrogate
model of the $0.95$-quantile of the outcome $Y$ of the experiments with both
technical systems.

As discussed in the last section in case of the MAFDS, the resulting improved
surrogate models are (almost) equal for both computer models and thus the
resulting quantile estimate is (almost equal).

##  5 Comparison of computer models

In this section methods which can be used to compare different computer
simulations are described.

###  5.1 Comparison of computer models via cumulative distribution functions

Roy and Oberkampf, (2011) introduced a method to detect model form
uncertainty. The term model form is in their context the result of all
assumptions, conceptualizations, abstractions, approximations, and
mathematical formulations on which the model relies. Consequently an imperfect
model is from their perspective the result of wrong assumptions in the model
form. To detect uncertainty in the model form, they proposed the so-called
Area Validation Metric (AVM). The metric compares the cumulative density
function of numerical simulated data and of experimental measurements. It uses
a sample of experimental and computer simulated data

$Y_{1},\ldots,Y_{n}\quad\text{and}\quad m(X_{n+1}),\ldots,m(X_{n+L_{n}})$ (19)

of size $n\in\mathbb{N}$ and $L_{n}\in\mathbb{N}$, where as described in
Section 3, $m\colon\mathbb{R}^{d}\to$ is a computer simulation which is
designed to emulate an experiment with a technical system and
$X\in\mathbb{R}^{d}$ is the corresponding input for the experiment and
$Y\in\mathbb{R}$ is the output. Based on (19) the cumulative density functions
are estimated by the empirical distribution function defined by

$\hat{F}_{Y}(t)=\frac{1}{n}\sum_{i=1}^{n}I_{(-\infty,Y_{i}]}(t)\quad\text{and}\quad\hat{F}_{m(X)}(t)=\frac{1}{L_{n}}\sum_{i=n+1}^{n+L_{n}}I_{(-\infty,m(X_{i})]}(t).$
(20)

Finally the AVM is calculated by

$\int|\hat{F}_{Y}(t)-\hat{F}_{m(X)}(t)|dt.$ (21)

###  5.2 A Bayesian approach for quantifying the computer model error

Kennedy and O’Hagan, (2001) proposed a Bayesian method for the calibration of
computer models in case that the underlying mathematical model is imperfect.
They assume that their data is given by

$Y=m(\theta,x)+\delta(x)+\epsilon$ (22)

where $m(\cdot,\cdot)$ describes the computer model depending on some
parameter $\theta$ and some vector $x$ describing uncertain factors of the
reality, where $\delta$ is the model error and $\epsilon$ is a normally
distributed error with expectation zero which contains for instance
measurement errors. They propose Bayesian techniques for the choice of the
optimal parameter value $\theta^{*}$ (so called calibration), which also model
the model error $\delta$ and are therefore able to quantify the model error of
the computer model. In the sequel we assume that our system is already
calibrated in the sense that our computer model is based on some physical
model for which we have adjusted the parameters already to reality based on
separate experiments with our technical system. In this case the model of
Kennedy and O’Hagan, (2001) simplifies to

$Y=m(x)+\delta(x)+\epsilon.$ (23)

Here $m(\cdot)$ is the (already calibrated) computer model, $\delta$ is the
discrepancy term of this model, and $\epsilon$ is a random normally
distributed error with expectation zero and variance $\lambda$. In addition
$x$ is (a possible random) value from $\mathbb{R}^{d}$ describing input
parameters of the experiment with the technical system (e.g., drop height in
experiments with the MAFDS, cf. Subsection 2.1).

From model (23) it is assumed that $n$ observations are observed, i.e., the
observed data is

$(X_{1},m(X_{1}),Y_{1}),\dots,(X_{n},m(X_{n}),Y_{n})$ (24)

where

$Y_{i}=m(X_{i})+\delta(X_{i})+\epsilon_{i}\quad(i=1,\dots,n).$ (25)

Here $\epsilon,\epsilon_{1},\dots,\epsilon_{n}$ are independent and
identically distributed. Kennedy and O’Hagan, (2001) uses a Gaussian process
with mean function $\mu(\cdot):\mathbb{R}^{d}\rightarrow\mathbb{R}$ and
covariance function
$c(\cdot,\cdot):\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}_{+}$
to model the discrepancy term $\delta(\cdot)$. This means that for any
$k\in\mathbb{N}$ and any $X_{1},\dots,X_{k}\in\mathbb{R}^{d}$ the vector

$(\delta(X_{1}),\delta(X_{2}),\dots,\delta(X_{k}))^{T}$

is multivariate normally distributed with mean vector
$(\mu(X_{1}),\dots,\mu(X_{k}))^{T}$ and covariance matrix
$(c(X_{i},X_{j}))_{1\leq i,j\leq k}$. The random error is chosen normally
distributed with mean zero and variance $\lambda$. So as soon as $\lambda$,
the mean function $\mu(\cdot)$ and the covariance function $c(\cdot,\cdot)$
are chosen, the model (23) is completely specified for given $x$.

In order to specify $\mu(\cdot)$ and $c(\cdot,\cdot)$, Kennedy and O’Hagan,
(2001) propose to use a hierarchical Bayesian model. Here the choice of
$\mu(\cdot)$ and $c(\cdot,\cdot)$ could be based on prior knowledge about the
technical system, but since this is difficult to incorporate in this model
they suggest to use the following simple choice: First prior densities for the
parameters of the hierarchical Bayesian model are chosen from a parametric
model. Therefore prior knowledge (of the technical system) can be used to
choose these densities, or if such knowledge is not available, normal
distributions or a choice of the form

$p(t)=c\cdot\frac{1}{t},$ (26)

with parameter $c\in\mathbb{R}_{+}$ and where $p(t)$ must be set to zero for
large $t$, is recommended.

In a second step they choose $\mu(\cdot)$ as a constant function

$\mu(x)=\beta\quad(x\in\mathbb{R}^{d})$ (27)

and use for
$c(\cdot,\cdot):\mathbb{R}^{d}\times\mathbb{R}^{d}\rightarrow\mathbb{R}$ the
definition $c(z_{1},z_{2})=\sigma^{2}R(z_{1}.z_{2})$, where

$R(z_{1},z_{2})=\exp\left(-\sum_{j=1}^{d}\omega_{j}\cdot(z_{1}^{(j)}-z_{2}^{(j)})^{2}\right)\quad(z_{1},z_{2}\in\mathbb{R}^{d}).$
(28)

As prior distributions they recommend normal distributions for $\lambda$ and
$\beta$, with mean and variance parameters
$\mu_{\lambda},\sigma^{2}_{\lambda},\mu_{\beta}$ and $\sigma^{2}_{\beta}$, and
for all other parameters they recommend to use (26), with parameters
$c_{\sigma^{2}},c_{\omega_{1}},\ldots,c_{\omega_{d}}$. Thus the hierarchical
Bayesian model uses parameters

${\cal P}^{(1)}=\\{\lambda,\beta,\sigma,\omega_{1},\ldots,\omega_{d}\\},$ (29)

where the vector of parameters is a realization of a random variable

$(\Lambda,B,\Sigma,\Omega_{1},\ldots,\Omega_{d}),$ (30)

with parametric densities

$p_{\Lambda}(\cdot),\ldots,p_{\Omega_{d}}(\cdot),$ (31)

with distribution parameters

${\cal
P}^{(2)}=\\{\mu_{\Lambda},\sigma^{2}_{\Lambda},\mu_{B},\sigma^{2}_{B},c_{\Sigma^{2}},c_{\Omega_{1}},\ldots,c_{\Omega_{d}}\\}.$
(32)

Kennedy and O’Hagan, (2001) recommended to choose fixed values for the
parameters in ${\cal P}^{(2)}$ whenever possible, e.g. from information about
the physical system, to reduce complexity.

The next step is to compute the posterior density of
$(\Lambda,B,\Sigma,\Omega_{1},\ldots,\Omega_{d})$ given the data (24). From
the above definitions we get that the conditional density of
$(Y_{1},\dots,Y_{n})$ given
$(\Lambda,B,\Sigma,\Omega_{1},\ldots,\Omega_{d})=(\lambda,\beta,\sigma,\omega_{1},\ldots,\omega_{d})$
is given by

$\displaystyle
f_{(Y_{1},\dots,Y_{n})|(\Lambda,B,\Sigma,\Omega_{1},\ldots,\Omega_{d})=(\lambda,\beta,\sigma,\omega_{1},\ldots,\omega_{d})}(y_{1},\dots,y_{n})=$
$\displaystyle\quad\frac{1}{\sqrt{(2\pi)^{n}\cdot
det(\Theta)}}\cdot\exp\Big{(}-\frac{1}{2}\cdot(y_{1}-m(X_{1})-\beta,\dots,y_{n}-m(X_{n})-\beta)\cdot\Theta^{-1}\cdot$
$\displaystyle\hskip
156.49014pt(y_{1}-m(X_{1})-\beta,\dots,y_{n}-m(X_{n})-\beta)^{T}\Big{)}$ (33)

where

$\Theta=\left(c(X_{k},X_{l})+\lambda\cdot 1_{\\{k=l\\}}\right)_{1\leq k,l\leq
n},$

from which we can conclude that we have

$\displaystyle f_{(Y_{1},\dots,Y_{n})}(y_{1},\dots,y_{n})=$ (34)
$\displaystyle\int\\!\\!\ldots\\!\\!\int
f_{(Y_{1},\dots,Y_{n})|(\Lambda,B,\Sigma,\Omega_{1},\ldots,\Omega_{d})=(\lambda,\beta,\sigma,\omega_{1},\ldots,\omega_{d})}(y_{1},\dots,y_{n})\cdot
p_{\Lambda}(\lambda)\cdots p_{\Omega_{d}}(\omega_{d})\,d\lambda\dots
d\omega_{d}$

and

$\displaystyle
f_{(\Lambda,B,\Sigma,\Omega_{1},\ldots,\Omega_{d})|(Y_{1},\dots,Y_{n})=(y_{1},\dots,y_{n})}(\lambda,\dots,\omega_{d})=$
$\displaystyle\hskip
14.22636pt\frac{f_{(Y_{1},\dots,Y_{n})|(\Lambda,B,\Sigma,\Omega_{1},\ldots,\Omega_{d})=(\lambda,\beta,\sigma,\omega_{1},\ldots,\omega_{d})}(y_{1},\dots,y_{n})\cdot
p_{\Lambda}(\lambda)\cdot\dots\cdot
p_{\Omega_{d}}(\omega_{d})}{f_{(Y_{1},\dots,Y_{n})}(y_{1},\dots,y_{n})},$ (35)

where $p_{\Lambda}(\cdot),\ldots,p_{\Omega_{d}}(\cdot)$, are the above
specified prior densities with the corresponding parameters ${\cal P}^{(2)}$ .
Note that the value of (5.2) also depends on the values of the distribution
parameters ${\cal P}^{(1)}$.

Kennedy and O’Hagan, (2001) propose next to estimate the parameters ${\cal
P}^{(1)}\cup{\cal P}^{(2)}$ by a maximum likelihood estimate and threat the
estimated values as if they were the real values, i.e. they proposed to
compute

$\displaystyle\hat{{\cal P}}^{(1)}\cup\hat{{\cal P}}^{(2)}=\arg\max_{{\cal
P}^{(1)}\cup{\cal
P}^{(2)}}f_{(\Lambda,B,\Sigma,\Omega_{1},\ldots,\Omega_{d})|(Y_{1},\dots,Y_{n})=(y_{1},\dots,y_{n})}(\lambda,\dots,\omega_{d}),$
(36)

where $\hat{{\cal P}}^{(1)}\cup\hat{{\cal P}}^{(2)}$ are the resulting
estimates of all parameters included in ${\cal P}^{(1)}\cup{\cal P}^{(2)}$. To
compute this estimate one can use a numerical (approximative) maximation of
the posterior density (where a maximum with respect to $2d+8$ variables has to
be computed). To reduce the complexity in solving (36) one can also complete
the square for $B$ in (5.2) to find

$\hat{\beta}=\frac{(1,\ldots,1)\cdot\Theta^{-1}\cdot(y_{1}-m(X_{1}),\ldots,y_{n}-m(X_{n}))^{T}}{(1,\ldots,1)\cdot\Theta^{-1}\cdot(1,\ldots,1)^{T}}.$

Then the above maximum has only be computed with respect to $2d+5$ variables,
since if $\hat{\beta}$ is known, we do not need to estimate it’s distribution
parameters $\mu_{B},\sigma_{B}^{2}$. Alternatively one can simplify the
computation by choosing for some of the above variables plausible estimates,
e.g., $B$ can be estimated by

$\hat{\beta}=\frac{1}{n}\sum_{i=1}^{n}(Y_{i}-m(X_{i})).$ (37)

As soon as we have these estimate we know the distribution of the error

$(Y_{i}-m(X_{i}))_{i=1,\dots,n}=(\delta(X_{i})+\epsilon_{i})_{i=1,\dots,n},$

which is a multivariate normal distribution with mean vector
$(\hat{\beta},\dots,\hat{\beta})^{T}$ and covariance matrix

$\left(\hat{\sigma}^{2}\cdot\exp\left(-\sum_{j=1}^{d}\hat{\omega}_{j}\cdot(X_{k}^{(j)}-X_{l}^{(j)})^{2}\right)+\hat{\lambda}\cdot
1_{\\{k=l\\}}\right)_{1\leq k,l\leq n}.$

From this result we can, e.g., estimate quantiles of the absolute model error

$(|\delta(X_{i})+\epsilon_{i}|)_{i=1,\dots,n}$

occurring in the data (24) via Monte Carlo in order to quantify the model
error of our simulation model.

###  5.3 Quantifying the computer model error via bootstrap

Wong et al., (2017) proposed the following method for quantifying the the
model error using a bootstrap approach:

First estimate a surrogate model $\hat{m}_{L_{n}}$ by (7) as described in
Subsection 4.1. Next to quantify the model error, $B$ bootstrap samples are
generated. Therefore compute the residuals as in Subsection 4.1 by (8) and
generate the bootstrap samples $(X_{1,b},\epsilon_{1,b})$, $\dots$,
$(X_{n,b},\epsilon_{n,b})$ $(b=1\ldots,B)$ by sampling randomly and with
replacement $n$ values of $(X_{1},\epsilon_{1})$, $\ldots$,
$(X_{n},\epsilon_{n})$.

Next we split the bootstrap samples in two parts, a learning and an evaluation
sample. Therefore we choose a $n_{l}\in\\{1,\ldots,n-1\\}$. Then we construct
$B$ estimators of the residuals $\hat{m}_{n_{l},b}^{\epsilon}$ $(b=1\ldots,B)$
by (9), where we replace $(X_{1},\epsilon_{1})$, $\ldots$,
$(X_{n},\epsilon_{n})$ by $(X_{1,b},\epsilon_{1,b})$, $\dots$,
$(X_{n_{l},b},\epsilon_{n_{l},b})$. Since we only use a part of the
experimental data to construct the residual estimators this can lead to
unsatisfying results, if the sample size is too small as discussed in
Subsection 4.1. In this case we recommend to use the weighted version of the
estimator instead, as defined in (11).

For each bootstrap sample we estimate the $\alpha$-quantile
$q_{|\epsilon|,\alpha}$ of the absolute model error on the remaining input
values evaluated with the bootstrap residual estimator. The $\alpha$-quantile
estimator is then defined by

$\hat{q}_{|\hat{m}_{n_{l},b}^{\epsilon}(X)|,n-n_{l},\alpha}=\min\left\\{y\in\mathbb{R}:\frac{1}{n-n_{l}}\sum_{i=n_{l}+1}^{n}I_{\\{|\hat{m}_{n_{l},b}^{\epsilon}(X_{i})|\leq
y\\}}\geq\alpha\right\\}.$ (38)

With this estimator one can compare two computer models by comparing the
median of their bootstrap $\alpha$-quantiles.

###  5.4 Application

In the following, we will illustrate and discuss the above described methods
by applying them to data of the MAFDS described in Subsection 2. Although it
might be possible to apply the methods to data of the piezo-beam, cf.
Subsection 2.2, we only have one computer model available for this system and
thus can not compare different computer models.

####  5.4.1 Comparison of computer models via cumulative distribution
functions

We apply the in Subsection 5.1 described method to experimental and simulated
output data from the MAFDS. Data generated by different computer models is
applied separately.

The integral in (21) is approximated via a Riemann sum calculated on
equidistant points. The resulting values for both systems can be found in
Table 3.

| MAFDS
---|---
computer model 1 | $0.00084$
computer model 2 | $0.00152$
Table 3:  Value of the area validation metric for both computer models of the
MAFDS.

A plot of the resulting estimates of the cumulative density functions for both
technical systems is shown in Figure 6.

Figure 6: Three estimates of the cumulative density function of the
experimental outcome of the MAFDS. An estimator based on experimental data
(black solid line), an estimator based on computer experiments with computer
model 1 (black dashed line) and computer model 2 (gray dashed line).

We can use the method to compare the two different computer models of the
MAFDS. The value of the AVM is smaller for computer model 1 and also the
estimated cumulative density function based on model 1 fits the estimate of
the cumulative density function based on the experimental data better, cf.
Figure 6. Thus based on the AVM one can conclude that model 1 predicts the
reality better than model 2.

The method can be used to compare different computer simulations, but
otherwise not much conclusions can be drawn, since the value of the AVM is
hard to interpret on its own.

####  5.4.2 A Bayesian approach for quantifying the computer model error

We apply the in Subsection 5.2 described method of Kennedy and O’Hagan, (2001)
on data of the MAFDS. The value of $\beta$ is estimated via (37). The
remaining parameters are estimated by the maximum likelihood estimator defined
by (36), where $\beta$ takes the fixed value estimated before and to
approximate the minimum we use a interior-point method, as implemented by the
MATLAB function fmincon.

We are interested in the $0.95$-quantiles of the $n=100$ absolute model
errors. Therefore we generate a sample of $10^{6}$ realizations of the $n=100$
model errors, using the estimated distribution parameters, calculate the
empirical $0.95$-quantiles and compute their median. The results are shown in
Table 4.

| MAFDS
---|---
computer model 1 | $0.161$
computer model 2 | $0.286$
Table 4:  Estimated $0.95$-quantiles of the $n=100$ maximal model errors for
both computer models of the MAFDS.

Figure 7 shows two box plots of $10^{6}$ estimates of the $0.95$-quantiles of
the $n=100$ absolute model errors for both technical systems.

Figure 7: Box plot without outliers of $10^{6}$ estimates for the
$0.95$-quantiles of the $n=100$ absolute model errors for both computer models
of the MAFDS.

In case of the MAFDS, we can use the method to compare the different computer
models and analyse their computer model error. The predicted computed model
error is smaller for computer model 1. Thus based on this method one can
conclude that model 1 predicts the reality better.

To conclude, the method can be used to compare different computer simulations
and estimate the model error.

####  5.4.3 Quantifying the computer model error via bootstrap

We apply bootstrap approach from Subsection 5.3 to data of the MAFDS and are
interested in the $0.95$-quantile of the absolute model error. We set $B=500$
and $n_{l}=10$ and define the estimate by (38). The resulting medians are
shown in Table 5.

| MAFDS
---|---
computer model 1 | $0.00138$
computer model 2 | $0.00326$
Table 5:  Median of the bootstrap estimates of the $0.95$-quantiles of the
$n=100$ maximal model errors for both computer models of the MAFDS.

Furthermore Figure 8 shows a box plot for both technical systems of the
$B=500$ bootstrap estimates of the $0.95$-quantiles.

Figure 8: Box plot without outliers for the $0.95$-quantiles of the $n=100$
absolute model errors for both computer models of the MAFDS. The box plot
contains $B=500$ bootstrap $0.95$-quantiles.

We can use the method to compare the different computer models and analyse
their computer model error. The predicted computed model error is smaller for
computer model 1. Thus based on this method one can conclude that model 1
predicts the reality better, with an estimated $0.95$-quantile of $0.00138$,

To conclude, the method can be used firstly to compare different computer
simulations and secondly to estimate the model error.

###  5.5 Comparison

All methods proved useful in the estimation of the computer model error,
although they differ in quality. The AVM of Roy and Oberkampf, (2011) is quite
simple. It measures the accuracy of computer models by estimates of their
cumulative distribution functions and compares these estimates to estimates of
the cumulative distribution function based on experimental data. The
difference is measured by the area between these estimates. But the result is
hard to interpret.

Using the method of Kennedy and O’Hagan, (2001) to quantify the model error,
it is necessary to compute a maximum likelihood estimate. This step is quite
unstable since one has to invert matrices, which leads in our examples to an
overestimation of the computer model error. In fact in case of the MAFDS the
maximal absolute model error of computer model 1 calculated on the $100$
experimental data is $0.0023$ and $0.0052$ for computer model 2, which is $70$
and $50$ times smaller than the estimated $0.95$-quantile of the model error,
cf. Table 4.

In the applications the method of Wong et al., (2017) proved the most
accurate. In fact, if we use the $100$ available data points to calculate the
empirical $0.95$-quantiles of the absolute model error for both computer
models of the MAFDS the resulting values are $0.0014$ for computer model 1 and
$0.0033$ for computer model 2. These estimates are quite close to the by the
method estimated values.

##  6 Quantifying the influence of the computer model error

In the following section, two different methods which analyse the influence of
the computer model error on density and quantile estimates are described. The
first method aims to estimate a confidence interval for quantiles. The second
method aims on estimating a confidence band of the density. With the aid of
these methods, it is also possible to compare the influence of the error of
computer models.

###  6.1 Confidence intervals for quantiles

Kohler et al., (2018) proposed a method to construct confidence intervals on
quantiles related to experiments with a technical system. The aim of the
method is to construct an estimator of the confidence interval of the
$\alpha$-quantile $q_{Y,\alpha}$ defined by

$q_{Y,\alpha}=\min\\{y\in\mathbb{R}\colon G_{Y}(y)\geq\alpha\\},$ (39)

where

$G_{Y}\colon\mathbb{R}\to[0,1],\quad G_{Y}(y)={\mathbf{P}}\\{Y\leq y\\},$ (40)

i.e. for a given $\delta\in(0,1)$ they want to construct a (random) interval
$\hat{C}_{n}$ such that

${\mathbf{P}}\\{q_{Y,\alpha}\in\hat{C}_{n}\\}\geq 1-\delta$ (41)

holds.

To do this, in a first step a surrogate model $\hat{m}_{L_{n}}$ of $m$ as in
Subsection 4.1 is defined by (7). Then to construct a confidence interval a
quantile estimator is constructed. Therefore define an estimator for the
cumulative density function of $\hat{m}_{L_{n}}(X)$ by

$\hat{G}_{\hat{m}_{L_{n}}(X),N_{n}}(y)=\frac{1}{N_{n}}\sum_{i=1}^{N_{n}}I_{(-\infty,\hat{m}_{L_{n}}(X_{n+L_{n}+i})]}(y)$
(42)

and define the corresponding plug in estimate by

$\hat{q}_{\hat{m}_{L_{n}}(X),N_{n},\alpha}=\min\\{y\in\mathbb{R}:\hat{G}_{\hat{m}_{L_{n}}(X),N_{n}}(y)\geq\alpha\\}$
(43)

for $\alpha\in(0,1)$. Next choose some suitable
$\varDelta\delta\in(0,\delta)$, set

$\displaystyle\hat{\beta}_{n}$ $\displaystyle=$
$\displaystyle\max_{i=1,\ldots,n}|Y_{i}-\hat{m}_{L_{n}}(X_{i})|,$
$\displaystyle N_{1}(n)$ $\displaystyle=$ $\displaystyle
N_{n}\cdot\left(\alpha-\sqrt{-\frac{\log(\varDelta\delta/2)}{2N_{n}}}\right),$
$\displaystyle N_{2}(n)$ $\displaystyle=$ $\displaystyle
N_{n}\cdot\left(\alpha+\sqrt{-\frac{\log(\varDelta\delta/2)}{2N_{n}}}\right),$

and choose $\epsilon_{n}\in(0,1)$ such that
$(1-\epsilon_{n})^{n}<(\delta-\varDelta\delta)$ holds and

$\epsilon_{n}+\gamma_{n}=\epsilon_{n}+\sqrt{-\frac{\log(\delta-\varDelta\delta-(1-\epsilon_{n})^{n})}{2N_{n}}}$
(44)

is minimal.

Finally the $\delta$ confidence interval for $q_{Y,\alpha}$ is defined by

$\hat{C}_{n}=\left[\hat{q}_{\hat{m}_{L_{n}}(X),N_{n},\frac{N_{1}(n)}{N_{n}}-\epsilon_{n}-\gamma_{n}}-\hat{\beta}_{n},\hat{q}_{\hat{m}_{L_{n}}(X),N_{n},\frac{N_{2}(n)}{N_{n}}+\epsilon_{n}+\gamma_{n}}+\hat{\beta}_{n}\right].$
(45)

###  6.2 Confidence bands for densities

Kohler and Krzyżak, (2019) proposed a method to construct confidence bands for
the density of $Y$. The aim of their method is to construct for a finite
sample size $n$ lower and upper bounds $\hat{g}_{n}^{(lower)}$ and
$\hat{g}_{n}^{(upper)}$ of $g$ satisfying

$\int_{I}\hat{g}_{n}^{(lower)}(x)\,dx\leq\int_{I}g(x)\,dx\leq\int_{I}\hat{g}_{n}^{(upper)}\,dx$
(46)

for all intervals $I$ with length $|I|\geq\kappa_{n}>0$, with probability
$1-\delta$, where $\delta,\kappa_{n}>0$ are given values. To define the
estimate the following notation is used: For $\beta>0$ define

$I^{\beta}=\\{x\in\mathbb{R}:[x-\beta,x+\beta]\cap I\neq\emptyset\\}$

and

$I_{\beta}=\\{x\in\mathbb{R}:[x-\beta,x+\beta]\subseteq I\\}.$

The empirical distribution $\hat{\mu}_{\hat{m}_{L_{n}(X),N_{n}}}$ of
$\hat{m}_{L_{n}}(X)$ is defined by

$\hat{\mu}_{\hat{m}_{L_{n}}(X),N_{n}}(A)=\frac{1}{N_{n}}\sum_{i=n+L_{n}+1}^{n+L_{n}+N_{n}}I_{A}(\hat{m}_{L_{n}}(X_{i}))\quad(A\subseteq\mathbb{R}).$

In a first step a surrogate model $\hat{m}_{L_{n}}$ of $m$ is defined as in
Subsection 4.1 by (7). Then construct a kernel density estimate based on a
sample of $\hat{m}_{L_{n}}(X)$, i.e. the estimator is defined by

$\hat{f}_{\hat{m}_{L_{n}}(X),N_{n},h_{N_{n}}}(y)=\frac{1}{N_{n}\cdot
h_{N_{n}}}\cdot\sum_{i=n+L_{n}+1}^{n+L_{n}+N_{n}}K\left(\frac{y-\hat{m}_{L_{n}}(X_{i})}{h_{N_{n}}}\right),$

for some bandwidth $h_{N_{n}}>0$ and some kernel
$K\colon\mathbb{R}\to\mathbb{R}$.

Set $\hat{\beta}=\max_{i=1,\ldots,n}|Y_{i}-\hat{m}_{L_{n}}(X_{i})|$ and choose
$\epsilon_{n}\in(0,1)$ such that $(1-\epsilon_{n})^{n}+2/N_{n}^{2}<\delta$
holds for a given $\delta\in(0,1)$ and that

$\gamma_{n}+\epsilon_{n}=\sqrt{\frac{-\ln(\delta-2/N_{n}^{2}-(1-\epsilon)^{n})}{2\cdot
N_{n}}}+\epsilon_{n}$ (47)

is minimal. Finally for some given $\kappa_{n}>0$ the upper and lower bounds
on the density $g$ are defined by

$\displaystyle\hat{g}_{n}^{(upper)}(y)$ $\displaystyle=$
$\displaystyle\hat{f}_{\hat{m}_{L_{n}}(X),N_{n},h_{N_{n}}}(y)+\frac{\epsilon_{n}+\gamma_{n}+\frac{2\cdot\sqrt{\log
N_{n}}}{\sqrt{N_{n}}}}{\kappa_{n}}$
$\displaystyle+\frac{1}{\kappa_{n}}\cdot\sup_{\begin{subarray}{c}J\text{
interval, }y\in J,\\\
|J|>\kappa_{n}\end{subarray}}\Bigg{(}\hat{\mu}_{\hat{m}_{L_{n}}(X),N_{n}}(J^{\hat{\beta}})-\int_{J}\hat{f}_{\hat{m}_{L_{n}}(X),N_{n},h_{N_{n}}}(t)dt\Bigg{)}$

and

$\displaystyle\hat{g}_{n}^{(lower)}(y)$ $\displaystyle=$
$\displaystyle\max\Bigg{\\{}\hat{f}_{\hat{m}_{L_{n}}(X),N_{n},h_{N_{n}}}(y)-\frac{\epsilon_{n}+\gamma_{n}+\frac{2\cdot\sqrt{\log
N_{n}}}{\sqrt{N_{n}}}}{\kappa_{n}}$
$\displaystyle-\frac{1}{\kappa_{n}}\cdot\sup_{\begin{subarray}{c}J\text{
interval, }y\in J,\\\
|J|>\kappa_{n}\end{subarray}}\Bigg{(}\int_{J}\hat{f}_{\hat{m}_{L_{n}}(X),N_{n},h_{N_{n}}}(t)dt-\hat{\mu}_{\hat{m}_{L_{n}}(X),N_{n}}(J_{\hat{\beta}})\Bigg{)},0\Bigg{\\}}.$

###  6.3 Application

In the following, we will illustrate and discuss the above described methods
by applying them to data of the technical systems described in Section 2.

####  6.3.1 Confidence intervals for quantiles

We apply the confidence interval estimator for quantiles from Subsection 6.1
to data of the MAFDS and are interested in the $95\%$ (i.e. $\delta=0.05$)
confidence interval of the $0.95$-quantile. The minimization problem in (44)
is approximated by a interior-point method, as implemented by the MATLAB
function fmincon.

The resulting confidence intervals for the two different computer models are
shown in Table 6.

| MAFDS
---|---
computer model 1 | $[0.0826,0.0890]$
computer model 2 | $[0.0786,0.0900]$
Table 6:  $95\%$ confidence intervals of the $0.95$-quantile of both computer
models of the MAFDS.

In case of the MAFDS, we can use the method to estimate confidence intervals
as above which yield more information about the experiment with the technical
system. Furthermore it is possible to compare different computer models by
analyzing the influence of the computer model error on the estimated
confidence intervals. The predicted confidence interval is narrower for
computer model 1 and since we used the same experimental data for both models,
the computer model error must be larger for computer model 2.

In case of the piezo-beam, there are only $10$ experimental data values
available and thus it is not possible to apply the above method for a
reasonable $\delta$ and $\alpha$.

To conclude, the method can be used to estimate a confidence interval of some
$\alpha$-quantile. It is only possible to apply it in case a sufficiently
large sample is available. In fact for the chosen values of $\alpha$ and
$\delta$, $n$ needs to be greater than $60$. For the method to work with $10$
data points and $\alpha=0.95$, one would need to increase $\delta$ to $0.63$.
Furthermore, the method can be used to analyse the influence of the computer
model error on the estimation of the confidence interval. Since the model
error directly influences the width of the confidence interval, one can also
compare different computer models.

####  6.3.2 Confidence bands for densities

We apply the confidence band estimator from Subsection 6.2 to data of both
technical systems. The minimization problem in (47) is approximated by a
interior-point method, as implemented by the MATLAB function fmincon. In both
cases we restrict the method to a finite interval. In case of the MAFDS we
only consider the interval $[0.075,0.09]$ and in case if the piezo-beam we
consider the interval $[12,15]$. We choose $\kappa_{n}$ as small as possible
to get a reasonable confidence band. In case of the MAFDS this means we set
$\kappa_{n}=0.005$ and in case of the piezo-beam we set $\kappa_{n}=2$. The
value of $\delta$ is set to $0.05$. Kohler and Krzyżak, (2019) pointed out
that the confidence band holds simultaneously for all bandwidths with
probability of at least $1-\delta$. I.e. we can calculate the confidence band
for a multitude of bandwidths and use pointwise the minimal value for the
upper band and the maximal value for the lower band. For computer model 1 of
the MAFDS we use $h_{N_{n}}\in\\{0.00005,0.0001,0.0005,0.0025,0.004,0.006\\}$,
for computer model 2 we use
$h_{N_{n}}\in\\{0.00001,0.0001,0.0005,0.0008,0.001\\}$ and for the piezo-beam
we use $h_{N_{n}}\in\\{0.005,0.05,0.1,0.2,0.3\\}$. As kernel function
$K\colon\mathbb{R}\to\mathbb{R}$ we use the naive kernel
$K(x)=\frac{1}{2}\cdot I_{[-1,1]}(x)$.

The results are displayed in Figure 9.

Figure 9:  Upper plots are corresponding to computer model 1 and 2 of the
MAFDS, cf. Subsection 2.1. Lower plot to piezo-beam, cf. Subsection 2.2. Each
plot contains density estimator (dotted line) based on a surrogate model, cf.
Subsection 4.2 and the upper and lower confidence bands
$\hat{g}_{n}^{(upper)}$ and $\hat{g}_{n}^{(lower)}$ (upper and lower solid
lines).

In case of the MAFDS, we can use the method to estimate a confidence band of
the density of $Y$ which yields more information about the experiment with the
technical system. Furthermore we can compare different computer models and
analyse their computer model error. Here the predicted confidence band is
narrower for computer model 1, hence the error of computer model 1 has a
smaller influence on the estimation of the density than the error of computer
model 2.

In case of the piezo-beam, we can use the method to estimate a confidence band
of the density of $Y$. This information can be used in an analysis or the next
design stage of the technical system.

To conclude, the method can be used to estimate a confidence band of the
density of the outcome of an experiment with a technical system. Furthermore,
the method can be used to analyse the influence of the computer model error on
the estimation of the confidence band. Since the model error influences the
width of the confidence interval, one can also compare different computer
models. Note that to achieve a reasonable confidence band in case of the
piezo-beam, it was necessary to set $\kappa_{n}$ to a rather large value, so
the results here are hard to interpret.

###  6.4 Comparison

The first method can be used to estimate confidence intervals for quantiles,
the second to estimate confidence bands for densities. By estimating a
confidence band for a density one gains more information about the technical
system. The problem is, that no real confidence band in a sup norm sense is
estimated, i.e. as stated by (46), the estimated confidence band only contains
the real density on intervals which are longer than $\kappa_{n}$ with
probability $1-\delta$. Thus it is quite hard to interpret the result.

On the other side, the estimated confidence intervals of quantiles do not
yield a lot information about the experimental outcome, but it is simpler to
analyse the influence of the model error and compare different computer models
as illustrated above.

##  7 Acknowledgment

Funded by the Deutsche Forschungsgemeinschaft (DFG, German Research
Foundation) - Projektnummer 57157498 - SFB 805.

## References

  * Bayarri et al., (2007) Bayarri, M. J., Berger, J. O., Cafeo, J., Garcia-Donato, G., Liu, F., Palomo, J., Parthasarathy, R. J., Paulo, R., Sacks, J., and Walsh, D. (2007).  Computer model validation with functional output.  Ann. Statist., 35(5):1874–1906.
  * Bichon et al., (2008) Bichon, B. J., Eldred, M. S., Swiler, L. P., Mahadevan, S., and McFarland, J. M. (2008).  Efficient global reliability analysis for nonlinear implicit performance functions.  AIAA Journal, 46(10):2459–2468.
  * Bott et al., (2015) Bott, A.-K., Felber, T., and Kohler, M. (2015).  Estimation of a density in a simulation model.  J. Nonparametr. Stat., 27(3):271–285.
  * Bourinet et al., (2011) Bourinet, J.-M., Deheeger, F., and Lemaire, M. (2011).  Assessing small failure probabilities by combined subset simulation and support vector machines.  Structural Safety, 33(6):343 – 353.
  * Bucher and Bourgund, (1990) Bucher, C. and Bourgund, U. (1990).  A fast and efficient response surface approach for structural reliability problems.  Structural Safety, 7(1):57 – 66.
  * Damblin et al., (2018) Damblin, G., Barbillon, P., Keller, M., Pasanisi, A., and Parent, E. (2018).  Adaptive numerical designs for the calibration of computer codes.  SIAM/ASA J. Uncertain. Quantif., 6(1):151–179.
  * Das and Zheng, (2000) Das, P. and Zheng, Y. (2000).  Cumulative formation of response surface and its use in reliability analysis.  Probabilistic Engineering Mechanics, 15(4):309 – 315.
  * Deheeger and Lemaire, (2010) Deheeger, F. and Lemaire, M. (2010).  Support vector machine for efficient subset simulations: 2smart method.  In Proceedings of the 10th Internatinal Conference on Applications of Statistics and Probability in Civil Engineering (ICASP10). Tokyo, Japan.
  * Devroye, (1986) Devroye, L. (1986).  Non-uniform random variate generation.  Springer-Verlag, New York.
  * Devroye et al., (2013) Devroye, L., Felber, T., and Kohler, M. (2013).  Estimation of a density using real and artificial data.  IEEE Trans. Inform. Theory, 59(3):1917–1928.
  * Enss et al., (2016) Enss, G. C., Kohler, M., Krzyżak, A., and Platz, R. (2016).  Nonparametric quantile estimation based on surrogate models.  IEEE Trans. Inform. Theory, 62(10):5727–5739.
  * Fang et al., (2006) Fang, K.-T., Li, R., and Sudjianto, A. (2006).  Design and modeling for computer experiments.  Chapman & Hall/CRC Computer Science and Data Analysis Series. Chapman & Hall/CRC, Boca Raton, FL.
  * (13) Felber, T., Kohler, M., and Krzyżak, A. (2015a).  Adaptive density estimation based on real and artificial data.  J. Nonparametr. Stat., 27(1):1–18.
  * (14) Felber, T., Kohler, M., and Krzyżak, A. (2015b).  Adaptive density estimation from data with small measurement errors.  IEEE Trans. Inform. Theory, 61(6):3446–3456.
  * Feldmann and Platz, (2019) Feldmann, R. and Platz, R. (2019).  Assessing fodel form uncertainty for a suspension strut using gaussian processes.  In UNCECOMP 2019; 3rd International Conference on Uncertainty Quantification in Computational Sciences and Engineering, Crete, Greece.
  * Götz, (2019) Götz, B. (2019).  Evaluation of uncertainty in the vibration attenuation with shunted piezoelectric transducers integrated in a beam-column support.  PhD thesis, Technische Universität.
  * Götz et al., (2018) Götz, B., Kersting, S., and Kohler, M. (2018).  Estimation of an improved surrogate model in uncertainty quantification by neural networks.  Submitted for publication.
  * Götz et al., (2016) Götz, B., Schaeffner, M., Platz, R., and Melz, T. (2016).  Lateral vibration attenuation of a beam with circular cross-section by a support with integrated piezoelectric transducers shunted to negative capacitances.  Smart Materials and Structures, 25(9):095045.
  * Gu and Wang, (2018) Gu, M. and Wang, L. (2018).  Scaled Gaussian stochastic process for computer model calibration and prediction.  SIAM/ASA J. Uncertain. Quantif., 6(4):1555–1583.
  * Han et al., (2009) Han, G., Santner, T. J., and Rawlinson, J. J. (2009).  Simultaneous determination of tuning and calibration parameters for computer experiments.  Technometrics, 51(4):464–474.
  * Higdon et al., (2013) Higdon, D., Gattiker, J., Lawrence, E., Jackson, C., Tobis, M., Pratola, M., Habib, S., Heitmann, K., and Price, S. (2013).  Computer model calibration using the ensemble Kalman filter.  Technometrics, 55(4):488–500.
  * Hurtado, (2004) Hurtado, J. E. (2004).  Structural reliability : statistical learning perspectives, volume Vol. 17 of Lecture notes in applied and computational mechanics.  Berlin [u.a.].
  * Kaymaz, (2005) Kaymaz, I. (2005).  Application of kriging method to structural reliability problems.  Structural Safety, 27(2):133–151.
  * Kennedy and O’Hagan, (2001) Kennedy, M. C. and O’Hagan, A. (2001).  Bayesian calibration of computer models (with discussion).  J. R. Stat. Soc. Ser. B Stat. Methodol., 63(3):425–464.
  * Kersting and Kohler, (2019) Kersting, S. and Kohler, M. (2019).  Uncertainty quantification based on (imperfect) simulation models with estimated input distributions.  Submitted for publication.
  * Kim and Na, (1997) Kim, S.-H. and Na, S.-W. (1997).  Response surface method using vector projected sampling points.  Structural Safety, 19(1):3 – 19.  Asian-Pacific Symposium on Structural Reliability and Its Applications.
  * Kohler and Krzyżak, (2017) Kohler, M. and Krzyżak, A. (2017).  Improving a surrogate model in uncertainty quantification by real data.  Submitted for publication.
  * Kohler and Krzyżak, (2018) Kohler, M. and Krzyżak, A. (2018).  Adaptive estimation of quantiles in a simulation model.  IEEE Trans. Inform. Theory, 64(1):501–512.
  * Kohler and Krzyżak, (2019) Kohler, M. and Krzyżak, A. (2019).  Estimation of extreme quantiles in a simulation model.  J. Nonparametr. Stat., 31(2):393–419.
  * Kohler et al., (2018) Kohler, M., Krzyżak, A., Mallapur, S., and Platz, R. (2018).  Uncertainty quantification in case of imperfect models: a non-Bayesian approach.  Scand. J. Stat., 45(3):729–752.
  * Li et al., (2017) Li, S., Götz, B., Schaeffner, M., and Platz, R. (2017).  Approach to prove the efficiency of the monte carlo method combined with the elementary effect method to quantify uncertainty of a beam structure with piezo–elastic supports.  pages 441–455.
  * Mallapur and Platz, (2019) Mallapur, S. and Platz, R. (2019).  Uncertainty quantification in the mathematical modelling of a suspension strut using bayesian inference.  Mechanical Systems and Signal Processing, 118:158–170.
  * McKay et al., (1979) McKay, M. D., Beckman, R. J., and Conover, W. J. (1979).  A comparison of three methods for selecting values of input variables in the analysis of output from a computer code.  Technometrics, 21(2):239–245.
  * Papadrakakis and Lagaros, (2002) Papadrakakis, M. and Lagaros, N. D. (2002).  Reliability-based structural optimization using neural networks and monte carlo simulation.  Computer Methods in Applied Mechanics and Engineering, 191(32):3491 – 3507.
  * Parzen, (1962) Parzen, E. (1962).  On estimation of a probability density function and mode.  Ann. Math. Statist., 33:1065–1076.
  * Plumlee, (2017) Plumlee, M. (2017).  Bayesian calibration of inexact computer models.  J. Amer. Statist. Assoc., 112(519):1274–1285.
  * Plumlee, (2019) Plumlee, M. (2019).  Computer model calibration with confidence and consistency.  J. R. Stat. Soc. Ser. B. Stat. Methodol., 81(3):519–545.
  * Rosenblatt, (1956) Rosenblatt, M. (1956).  Remarks on some nonparametric estimates of a density function.  Ann. Math. Statist., 27:832–837.
  * Roy and Oberkampf, (2011) Roy, C. J. and Oberkampf, W. L. (2011).  A comprehensive framework for verification, validation, and uncertainty quantification in scientific computing.  Comput. Methods Appl. Mech. Engrg., 200(25-28):2131–2144.
  * Santner et al., (2018) Santner, T. J., Williams, B. J., and Notz, W. I. (2018).  The design and analysis of computer experiments.  Springer Series in Statistics. Springer, New York.  Second edition of [ MR2160708].
  * Tuo and Wu, (2015) Tuo, R. and Wu, C. F. J. (2015).  Efficient calibration for imperfect computer models.  Ann. Statist., 43(6):2331–2352.
  * Wang et al., (2009) Wang, S., Chen, W., and Tsui, K.-L. (2009).  Bayesian validation of computer models.  Technometrics, 51(4):439–451.
  * Wong et al., (2017) Wong, R. K. W., Storlie, C. B., and Lee, T. C. M. (2017).  A frequentist approach to computer model calibration.  J. R. Stat. Soc. Ser. B. Stat. Methodol., 79(2):635–648.

## Supplementary material

### Completion of the square for $\beta$

In the sequel we use the following notation

$\bar{y}=(y_{1},\ldots,y_{n}),\;\bar{m}=(m(x_{1}),\ldots,m(x_{n}))\text{ and
}\bar{1}=(1,\ldots,1),$

where $\bar{1}$ is also $n$-dimensional. In the conditional density

$\displaystyle
f_{(Y_{1},\dots,Y_{n})|(\Lambda,\dots,\Omega_{d})=(\lambda,\dots,\omega_{d})}(y_{1},\dots,y_{n})$
$\displaystyle=\frac{1}{\sqrt{(2\pi)^{n}\cdot
det(\Theta)}}\cdot\exp\left(-\frac{1}{2}\cdot(\bar{y}-\bar{m}-\beta\cdot\bar{1})\cdot\Theta^{-1}\cdot(\bar{y}-\bar{m}-\beta\cdot\bar{1})^{T}\right)$

we want to complete the square for $\beta$. Thus we want to complete the
square in the exponent. Since $\Theta$ is symmetric and positive definite we
have that $\Theta^{-1}$ is symmetric and positive definite as well. Hence

$\displaystyle\exp\left(-\frac{1}{2}\cdot(\bar{y}-\bar{m}-\beta\cdot\bar{1})\cdot\Theta^{-1}\cdot(\bar{y}-\bar{m}-\beta\cdot\bar{1})^{T}\right)$
$\displaystyle=\exp\left(-\frac{1}{2}\cdot\left(\beta^{2}\cdot\bar{1}\cdot\Theta^{-1}\cdot\bar{1}^{T}-2\cdot\beta\bar{1}\cdot\Theta^{-1}\cdot(\bar{y}-\bar{m})^{T}+(\bar{y}-\bar{m})\cdot\Theta^{-1}\cdot(\bar{y}-\bar{m})^{T}\right)\right)$
$\displaystyle=\exp\left(-\frac{1}{2}\cdot\left(\bar{1}\cdot\Theta^{-1}\cdot\bar{1}^{T}\left(\beta-\frac{\bar{1}\cdot\Theta^{-1}\cdot(\bar{y}-\bar{m})^{T}}{\bar{1}\cdot\Theta^{-1}\cdot\bar{1}^{T}}\right)^{2}-\frac{\bar{1}\cdot\Theta^{-1}\cdot(\bar{y}-\bar{m})^{T}}{\bar{1}\cdot\Theta^{-1}\cdot\bar{1}^{T}}\right)\right)\cdot$
$\displaystyle\qquad\exp\left(-\frac{1}{2}\cdot(\bar{y}-\bar{m})\cdot\Theta^{-1}\cdot(\bar{y}-\bar{m})^{T}\right)$

holds. Since $\Theta^{-1}$ is symmetric and positive definite

$\bar{1}\cdot\Theta^{-1}\cdot\bar{1}^{T}>0$

holds. Thus to maximize the exponent above for $\beta$

$\beta=\frac{\bar{1}\cdot\Theta^{-1}\cdot(\bar{y}-\bar{m})^{T}}{\bar{1}\cdot\Theta^{-1}\cdot\bar{1}^{T}}$

needs to hold.

### Density

Since the proposed density for the parameters

$p(t)=const\cdot\frac{1}{t}$

is no real density we propose to use

$p(t)=I_{[\epsilon,\exp(1/const)\cdot\epsilon]}\cdot const\cdot\frac{1}{t}$

instead.