Thomas Mühlenstädt

Laurent Carraro

Clément Chevalier

David Ginsbourger

Bertrand Iooss

Roland Jegou

Art B Owen

Clémentine Prieur

Marc Roelens

David Steinberg

Mom- Chil Ivanov

Hamad Ul Hassan

Leo Geppert

Max Wornowizki

André König

Uwe Ligges

Olaf Mersmann

Daniel Vogel

Alper Güner

Jörg Kolbe

Stefan Hess

Simon Neumärker

Malte Jastrow

Matthias Borowski

Thoralf Mildenberger

André Rehage

Nikolaus Rudak

Henrike Weinert

Ieva Zelo

Christine Exbrayat

Simon Ferber

Carina Hösl

Stefanie Krombacher

Natalina Külow

This work presents new sensitivity analysis methods, motivated by and framed around applications in the analysis of shape accuracy in sheet metal forming. Shape deformation after the removal of the forming punch is a serious problem in production processes, but can be simulated by finite elements methods (FEM).

. Furthermore, similar to the input variable screening by Sobol indices, we can use TIIs for interaction screening. This leads to a block-additive decomposition of the computer experiment. An inactive interaction indicates that the two input variables come from two groups of variables that have no interaction in common and thus are additive. This knowledge about a block-additive structure in the computer experiment can be applied in subsequent procedures. Metamodels like Kriging can be adapted to follow the block-additive structure and optimization can be simplified and parallelized when being performed in each block separately. We present several estimation methods for the TII and analyze their statistical properties. Methods for interaction screening, including a suitable thresholding algorithm for inactive interactions, are developed. In addition, derivative-based indices, which provide a faster-to-compute upper bound for the TII, are presented.

Usually the process parameters are kept on a constant level during the computer experiment. For the deep drawing analysis, finite element simulations have been developed that allow the user to change parameters during the forming process in order to further improve the accuracy of the analysis. Here new methods are necessary for the sensitivity analysis of those temporal changeable parameters. This work presents an idea

Introduction

This work deals with the topic of sensitivity analysis for computer experiments. In the situation that statistical experiments are too expensive, time-consuming, might harm lives, or are simply impossible to perform, computer simulations can sometimes be used as a replacement. These computer experiments are complex mathematical models constructed from known physical relations, e.g. in engineering or environmental problems, or chemical and medical experiments, where computer experiments are also known as in silico experiments. Specific examples for computer experiments are traffic simulation models [START_REF] Punzo | Sensitivity analysis of car-following models: methodology and application[END_REF], where the interest lies in observing those effects that lead to congestions, models of manufacturing processes like the deep drawing model which motivated this work, or highly complex weather models, used in the climate impact research [START_REF] Katz | Techniques for estimating uncertainty in climate change scenarios and impact studies[END_REF]. Computer experiments have become everyday procedures in those fields and are also rising in other areas like social sciences and economics. None the less, it has to be kept in mind that computer experiments are only approximations of the real process and should be validated by real experiments whenever possible.

From the statistical point of view, computer experiments have to be treated specially since they are usually deterministic, contain a high number of input variables and are highly complex in terms of nonlinear relations and interactions. The inherent mathematical equations are most often solved numerically which results in long-lasting 1 computations, so computer experiments can be very time-consuming. Methods have been developed for the design of computer experiments, above all space-filling designs like Latin hypercubes, maximin design, or entropy design [START_REF] Santner | The design and analysis of computer experiments[END_REF]. For deterministic modeling and prediction, advanced interpolation methods exist which respect the high complexity, with the most popular being the Gaussian process emulator, also called Kriging. Finally, there are also optimization procedures for computer experiments, which make use of the prediction models, one popular method being the Kriging-based EGO algorithm, which utilizes the fact that Kriging provides an estimation of the uncertainty of unobserved points.

In this work, the focus lies on a further broad field in the analysis of computer experiments, sensitivity analysis, which has importance in all steps beginning with the building of the computer experiment over its use and understanding to model building and model validation. It analyzes how sensitive the computer experiment responds to variations in the input.

The analysis of the response's sensitivity is of interest on its own, e.g. to answer the question what happens in a car crash simulation when the tire friction changes. Based on [START_REF] Saltelli | Wiley series in probability and statistics[END_REF], five further points of sensitivity analysis application can be named:

(1) To check, if a model resembles the real process, i.e. if sensitivities reflect expectations,

(2) to determine the most important input variables and optimal regions in the space of input variables for calibration studies,

(3) to select and rank input variables by importance, (4) to detect regions in the space of input variables for which the model variation is maximum, and

(5) to detect (groups of) input variables that interact with each other.

for an effective sensitivity analysis of functional input including design and graphical representation of the functional influences of the inputs. To keep the number of function evaluations as low as possible, a sequential algorithm is introduced that increases the accuracy of the functional sensitivity analysis with every step. The points (1), ( 2) and (3) can be handled by the methods for functional input variables.

The ideas of the sequential algorithm are further developed for the analysis of the support of scalar input variables, the support analysis. Regarding point (4) the aim is to analyze the sensitivity of the different regions of the variables' distribution support. This information gives closer insight into the impact of the variables in the process and can furthermore help in finding optimal support settings for screening, modeling, and sensitivity analysis. In this work, support index functions as well as visualization methods are suggested and their relations to other indices are derived.

The thesis is structured as follows. Chapter 2 gives an introduction to sensitivity analysis together with a panorama of the current state of research in the field with special focus on variance-based global sensitivity indices. Chapter 3 introduces the TII along with the complete methodology for sensitivity analysis of the interaction structure. Chapter 4 presents the approach of sensitivity analysis for functional input and in Chapter 5 the developed technique for support analysis is presented. All methods are then applied in deep drawing simulations in Chapter 6. The work is concluded with a summary and outlook in Chapter 7.

Chapter 2

Sensitivity Analysis

Background

In the field of computer experiments, sensitivity analysis generally explores the relationships between information flowing in and out of the experiment. More specifically, it studies how the uncertainty in the output of the computer experiment can be apportioned to different sources of uncertainty in the inputs with the direct aim of process understanding and input variable screening, but also calibration, metamodel building and finally optimization. An elaborate overview can be found in the standard work of sensitivity analysis by [START_REF] Saltelli | Wiley series in probability and statistics[END_REF] and in the comparison of sensitivity analysis techniques by [START_REF] Confalonieri | Comparison of sensitivity analysis techniques: A case study with the rice model warm[END_REF]. This chapter starts with a brief outline of the background of sensitivity analysis.

In the beginning, sensitivity analysis was based on standard statistical methods like scatter plots. There, the influence of a variable can be read from plots of the output against each input variable. Further methods were regression analysis, where regression coefficients qualify the linear sensitivity of each input (e.g. [START_REF] Kleijnen | Sensitivity analysis and related analyses: a review of some statistical techniques[END_REF]) and correlation analysis. A next step were one-at-a-time methods (e.g. [START_REF] Daniel | One-at-a-time plans[END_REF]).

Here, one input variable is changed while keeping all other at a constant level, thus information can be obtained only for that particular region in the input space. A similar approach came from chemical modeling (e.g. Chapter 5 of [START_REF] Saltelli | Wiley series in probability and statistics[END_REF]), which uses partial derivatives around a nominal value as local sensitivities. In contrast to those local methods, [START_REF] Schaibly | Study of the sensitivity of coupled reaction systems to uncertainties in rate coefficients. II Applications[END_REF] developed a global sensitivity analysis method called FAST (Fourier amplitude sensitivity test), which allows for all input variables to be varied at the same time. Here, the input variables are sampled in such a way that the amplitudes obtained by a Fourier analysis of the output can be interpreted as sensitivity indices of the input variables. Later on, Sobol ' (1993) introduced indices, which are now called Sobol indices, global sensitivity indices, which have been proved by [START_REF] Saltelli | An alternative way to compute Fourier amplitude sensitivity test (FAST)[END_REF] to predict the same quantities as the FAST indices. Because of their clear interpretation, unbiased estimation and their access to interactions, Sobol indices have been widely used and are being further developed continuously. Parallel to that, alternative sensitivity analysis methods were introduced like derivative-based indices [START_REF] Kucherenko | Monte Carlo evaluation of derivative-based global sensitivity measures[END_REF] and, based on oneat-a-time methods, very effective screening designs have been developed with the most important being Morris screening [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF].

The role of metamodels

The estimation of sensitivity indices usually requires a high number of model evaluations, as crude discrete integration methods are applied. Time-consuming computer experiments are therefore often replaced by a faster-to-evaluate second model, which is called a metamodel, surrogate model, emulator or response surface. Possible metamodels come from the large field of the modeling of computer experiments (e.g. [START_REF] Fang | Design and modeling for computer experiments[END_REF], a popular one being the Gaussian process model described below, but also polynomials, Polynomial Chaos Expansion, splines, or Neural Networks. As the metamodel is only an approximation of the real experiment, a metamodel error has to be taken into account and the accuracy of the model has to be assessed carefully, e.g. by leave-one-out cross validation.

Throughout this thesis, there will be no distinction between the true computer experiment and the metamodel. The term underlying model, symbolized by f , can either refer to the direct computer experiment or to a metamodel with negligible error. The underlying model f will be regarded as a black-box function, whose specific shape is inaccessible. However, there are some approaches that are able to utilize the specific properties of a metamodel. [START_REF] Blatman | Efficient computation of global sensitivity indices using sparse polynomial chaos expansions[END_REF] compute Sobol indices analytically from Sparse Polynomial Chaos Expansion models and [START_REF] Marrel | Calculations of Sobol indices for the Gaussian process metamodel[END_REF] define indices for the Gaussian process model that use the Gaussian process directly instead of the prediction function.

The Gaussian process model, also called Kriging as it was proposed by [START_REF] Krige | A statistical approach to some basic mine valuation problems on the witwatersrand[END_REF], is a standard tool in computer experiments for various reasons. It interpolates the data, which is important for deterministic computer codes. In addition, the Kriging model performs well in predicting, even for highly complex functions, and it gives a measure of uncertainty at unobserved points [START_REF] Fang | Design and modeling for computer experiments[END_REF]. The basic idea is to assume that the function is a realization of a Gaussian process with a linear trend as mean,

Y (x) = p =1 β f (x) + Z(x),
(2.1) where x = (x 1 , . . . , x d ) represents the vector of d input variables, f 1 (x), . . . , f p (x) are p known regression functions with β 1 , . . . , β p the corresponding parameters, and Z(.) is a Gaussian process with zero mean and covariance function, or kernel k. A particular class for the covariance kernel is the stationary family, which implies that k depends only on the difference between two different locations, k(x (1) -x (2) ) with k(.) = σ 2 R(.; θ), where σ 2 is the process variance, R the correlation function and θ a vector of covariance parameters.

As the Gaussian process model in (2.1) returns a Gaussian process instead of a function, the predictor has to be defined separately,

y(x * ) = p =1 β f (x * ) + k(x * ) K -1 (y -F β), (2.2)
with x * a new point at which to predict, K the covariance matrix at all data points, k(x * ) the covariance vector between the data and x * , and F the experimental matrix containing the trend values at all data points. Usually, the parameters β k , σ 2 and θ are estimated by numerical maximization of the likelihood and inserted into Equation (2.2).

In practice, the kernels are often modeled as tensor products of one-dimensional kernels,

k(h) = σ 2 d =1 g (h ; θ ), (2.3)
where popular one-dimensional covariance functions are the Gaussian,

g(h; θ) = exp - h 2 2θ 2 ,
and the Matérn 5/2,

g(h; θ) = 1 + √ 5|h| θ + 5h 2 3θ 2 exp - √ 5|h| θ .
An overview of possible covariance functions and their properties can be found in [START_REF] Rasmussen | Gaussian processes for machine learning[END_REF]. In the one-dimensional covariance functions, the parameter to be estimated, θ i , controls the covariance in the direction of input variable

x i . If all input variables are scaled equally, the estimates of θ i can be used as a first impression of the influence of the variables. A high covariance indicates a flat curve and thus a small influence of the variable whereas a low covariance indicates a higher fluctuation and a stronger influence.

Variance-based indices

The most popular global sensitivity analysis method, the so-called Sobol indices, is based on the variance of the function, decomposed into additive terms. One of the first to mention this decomposition was [START_REF] Hoeffding | A class of statistics with asymptotically normal distribution[END_REF], who used it to obtain independent random variables to study properties of U-statistics. [START_REF] Efron | The jackknife estimate of variance[END_REF] showed the uniqueness of the decomposition and Sobol ' (1993) revisited it in the context of sensitivity analysis. A generalization of the decomposition, the so-called high-dimensional model representation (HDMR), was introduced by [START_REF] Rabitz | Efficient input-output model representations[END_REF], who aimed at representing functions as a sum of low-dimensional components. Basing on so many (and more) contributors, the decomposition is known under a variety of names. In this work, it will be referred to as functional ANOVA (FANOVA) decomposition as it provides an ANOVA decomposition of the variance of the function.

Let us consider a black-box function Y = f (X), where X = (X 1 , . . . , X d ) is a vector

of independent random variables with distribution µ = µ 1 ⊗ • • • ⊗ µ d , and f is a d- dimensional function f : ∆ → R with f (X) ∈ L 2 (µ)
, where L 2 (µ) denotes the space of square-integrable functions with respect to the measure µ. The function can be decomposed into additive terms,

f (X) = f 0 + d i=1 f i (X i ) + i<j f i,j (X i , X j ) + • • • + f 1,...,d (X 1 , . . . , X d ).
(2.4)

The terms represent first-order effects (f i (X i )), second-order interactions (f i,j (X i , X j )), and all higher combinations of input variables. The decomposition is unique if all terms f I (X I ), I ⊂ {1, . . . , d} have zero mean, E(f I (X I )) = 0, I ⊆ {1, . . . , d} (2.5) and the conditional expectations fulfill the non-simplification conditions

E(f I (X I ) | X J ) = 0, J ⊂ I ⊆ (1, . . . , d) .
(2.6) From (2.5) and (2.6), it follows that the terms also have zero correlation,

E(f I (X I )f I (X I )) = 0, I = I .
The decomposition can be obtained by recursive integration,

f 0 = E(f (X)), f i (X i ) = E(f (X)|X i ) -f 0 , f i,j (X i , X j ) = E(f (X)|X i , X j ) -f i (X i ) -f j (X j ) -f 0 ,
and so on. By computing the variance of (2.4), an ANOVA-like variance decomposition is obtained, where each part quantifies the impact of the input variables on the response,

D = Var(f (X)) = Var(f 0 ) + d i=1 Var(f i (X i )) + i<j Var(f i,j (X i , X j )) + • • • + Var(f 1,...,d (X 1 , . . . , X d )).
(2.7)

Definition of various variance-based indices

The variances

D I = Var(f I (X I )) (2.8)
are known as unscaled Sobol indices [START_REF] Sobol | Sensitivity estimates for non linear mathematical models[END_REF]. The first-order Sobol index (for I ∈ {1, . . . , d}) is widely used as a sensitivity measure for quantifying the influence of firstorder effects. When I contains more than one input variable, the Sobol index quantifies the pure interaction influence of the variables indexed in I. As the estimation of the interactions of all possible variable combinations is usually laborious, other indices presented in the following have been developed for the assessment of interactions.

The Sobol index is often divided by the overall variance D, leading to the scaled Sobol index

S I = Var(f I (X I )) D = D I D .
(2.9) By this division by the overall variance, the index is normalized to fall between 0 and 1 and thus is easier to assess. The same applies for all indices introduced in the following.

For the sake of brevity, they will however be introduced in their unscaled versions since the scaling is equal for all methods considered.

An important extension of the Sobol index is the total sensitivity index D T I [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF]. Defined for a group of input variables X I for any I ⊆ {1, . . . , d}, it describes the influence of the variables including all interactions of any order that contain at least one of them. Thus, the influence on all orders instead of only the first is measured, which makes it more useful for the screening of input variables. The total sensitivity index is defined as the sum of all partial variances that contain at least one of the variables,

D T I = J∩I =∅ D J .
(2.10)

Another way to describe the influence of a group of input variables is the closed sensitivity index D C I , see e.g. [START_REF] Sobol | Sensitivity estimates for non linear mathematical models[END_REF],

D C I = J⊆I D J = Var (E[f (X)|X I ]) .
(2.11)

In contrast to total indices, interactions with variables not in X I are not included in this index, but all effects, first-order effects as well as interactions, caused by subsets of it. It is equal to the variance of the conditional expectation and is therefore also known under this name (or shorter VCE) in the literature [START_REF] Mckay | Evaluating prediction uncertainty[END_REF]. In its first-order version, it matches the Sobol index.

Sobol indices of higher order can be obtained as a sum of closed sensitivity indices of this order and lower,

D I = |I| M =1 (-1) |I|-M J⊆I, |J|=M D C J .
(2.12)

The total and the closed sensitivity index are related by (2.13) where -I denotes the complement of the set I. Through the formula of total variance, this leads to a further representation of the total sensitivity index,

D = D C -I + D T I ,
D T I = E (Var[f (X)|X -I ]) .

Monte Carlo estimation

As the computer experiments considered here are treated as black-box functions, the different indices cannot be computed directly, but have to be estimated. The most forward way to estimate the indices, that is the expected values resulting from (2.4), is by crude Monte Carlo integration. A high number n -typically not less than 1 000 × d -of random or quasi-random sampled Monte Carlo runs x (1) , . . . , x (n) is drawn from µ, the distribution of X, to approach the integral,

1 n n k=1 f (x (k) ) n→∞ ---→ f (X) dµ(X) = E(f (X)).
(2.14)

The size of the Monte Carlo sample is denoted by n throughout this work, except when its specific role is of interest. The approximation is unbiased and, following the strong law of large numbers, convergent with probability one [START_REF] Caflisch | Monte Carlo and quasi-Monte Carlo methods[END_REF].

To estimate the closed sensitivity indices, and via them through (2.12) the Sobol indices, we need a representation of the index that can directly be adapted to Monte Carlo integration. A common way, sometimes called the pick-and-freeze (or just pick-freeze) formula, is given by Sobol' (1993), (2.15) where here as well as in the following, Z stands for an independent copy of X. The constant term and overall variance can be obtained by

D C I = E [f (X)f (X I , Z -I )] -f 2 0 ,
f 0 = E(f (X)) and D = Var(f (X)).
The corresponding Monte Carlo estimate for two n × d-Monte Carlo samples x and z drawn from µ is straightforward.

pf D C I = 1 n n k=1 f (x (k) )f (x (k) I , z (k) -I ) -f 2 0 , (2.16) f 0 = 1 n n k=1 f (x (k) ),
(2.17)

D = Var (f (x)) . (2.18)
Alternative expressions, which also take the evaluations f (X I , Z -I ) into account in the estimation of f 0 and D and thus enable greater numerical stability, have been introduced by Monod et al. (2006, p. 86),

f 0 = E f (X) + f (X I , Z -I ) 2 , * f 0 = 1 n n k=1 f (x (k) ) + f (x (k) I , z (k) -I ) 2 , (2.19) and D = E f (X) 2 + f (X I , Z -I ) 2 2 -E f (X) + f (X I , Z -I ) 2 2 , * D = 1 n n k=1 f (x (k) ) 2 + f (x (k) I , z (k) -I ) 2 2 - * f 2 0 .
(2.20)

One problem of the pick-freeze estimator (2.16) is that its variance gets very large when f 2 0 is large in comparison to D C I . Sobol' (1993) therefore suggests to shift f by an amount close to f 0 . Another possibility is to avoid the subtraction of f 2 0 . In Owen (2013a), two such estimation strategies called Correlation 1 (Mauntz, 2002, Formula 18) and Correlation 2 are compared. They are based on the representations

D C I = E [f (X) (f (X I , Z -I ) -f (Z))] = E [(f (X) -f (Z I , X -I )) (f (X I , U -I ) -f (U ))] ,
with U a further independent copy of X. The corresponding estimates, using Monte Carlo samples x, z, and u from µ, are much more accurate for small closed sensitivity indices.

Cor1 D C I = 1 n n k=1 f (x (k) ) f (x (k) I , z (k) -I ) -f (z (k) ) , Cor2 D C I = 1 n n k=1 f (x (k) ) -f (z (k) I , x (k) -I ) f (x (k) I , u (k) -I ) -f (u (k) ) .
(2.21)

Though they require more -3 and 4 -vectors of function evaluations in the integral, no additional evaluations are necessary when applying the strategy of simultaneous estimation of all indices by [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF], described later in this Section.

D (overall variance) J D J Var(f (X)) E f (X) 2 -f 2 0 E f (X I ) 2 +f (X I ,Z -I ) 2 2 -f 2 0 D I (Sobol index) D I Var(f I (X I )) |I| M =1 (-1) |I|-M J⊆I, |J|=M D C I D C I (closed sensitivity index) J⊆I D J Var (E[f (X)|X I ]) E [f (X)f (X I , Z -I )] -f 2 0 , E [f (X) (f (X I , Z -I ) -f (Z))] , E [(f (X) -f (Z I , X -I )) (f (X I , U -I ) -f (U ))] D T I (total sensitivity index) J∩I =∅ D J E (Var[f (X)|X -I ]) D -D C I , 1 2 E (f (X) -f (X -I , Z I )) 2
Table 2.1: Overview of (unscaled) variance-based indices and their computations.

The total sensitivity index can on the one hand be obtained via the pick-freeze formula using the relation between closed and total sensitivity indices (2.13),

D T I = D -D C -I = D -E [f (X)f (Z I , X -I )] + f 2 0 . (2.22)
Another way, sometimes called the Jansen formula, which avoids the outer sum was first mentioned in Sobol ' (1993) and later improved in Jansen (1999),

D T I = 1 2 E (f (X) -f (Z I , X -I )) 2 .
(2.23)

The corresponding estimator

Jan D T I = 1 2n n k=1 f (x (k) ) -f (z (k) I , x (k) -I ) 2 .
is proved to be more efficient than (2.22) in Sobol' (2001, Theorem 4). A proof of its asymptotical efficiency will follow in Chapter 3.

An overview of all presented variance-based indices together with their estimation methods can be found in Tab. 2.1. (2002) presents a strategy to recycle runs when the full set of all d indices shall be estimated simultaneously, which is usually the case in applications. They exploit the fact that the two input data sets x and z can be used in more than one way in the estimation of an index. For instance the two pick-freeze formulas

Estimation of a full set of indices

Saltelli

1 n n k=1 f (x (k) )f (x (k) i,j , z (k)
-{i,j} ) -f 2 0 and

1 n n k=1 f (z (k) i , x (k) -i )f (x (k) j , z (k)
-j ) -f 2 0 are both estimators for D C i,j , as they both freeze the variables corresponding to {i, j}. Only the names are changed from x i and x j in the first formula to z i and x j in the second one. [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF] shows in his Theorem 1 that, using n(d + 2) evaluations

f (x), f (z), f (x 1 , z -1 ), . . . , f (x d , z -d ), it is possible to compute a) each first-order Sobol index D i = 1 n n k=1 f (x (k) )f (x (k) i , z (k) -i ) -f 2 0 , b) each total sensitivity index D T i = D -1 n n k=1 f (z (k) )f (x (k) i , z (k)
-i ) + f 2 0 , c) each (d -2)-order closed effect D C -{i,j} = 1 n n k=1 f (x (k) i , z (k)
-i )f (x (k) j , z (k)
-j ) -f 2 0 .

D and f 2 0 can be estimated from the evaluations as well, e.g. by (2.19) and (2.20). The method further allows for the estimation of all second-order total sensitivity indices and all the corresponding scaled indices.

If in addition n × d more evaluations f (z 1 , x -1 ), . . . , f (z d , x -d ) are available, it is possible to obtain second-order closed sensitivity indices as well as double estimates, which allow for more efficient estimation by taking the mean of both estimates. More precisely, we can compute a) an additional estimate of each first-order Sobol index

D i = 1 n n k=1 f (z (k) )f (z (k) i , x (k) -i ) -f 2 0 , b) an additional estimate of each total sensitivity index D T i = D -1 n n k=1 f (x (k) )f (z (k) i , x (k) -i ) + f 2 0 , c) an additional estimate of each (d -2)-order closed effect D C -{i,j} = 1 n n k=1 f (z (k) i , x (k) -i )f (z (k) j , x (k)
-j ) -f 2 0 , d) double estimates of each second-order closed effect

D C i,j = 1 n n k=1 f (x (k) i , z (k) -i )f (z (k) j , x (k) -j ) -f 2 0 and D * C i,j = 1 n n k=1 f (z (k) i , x (k) -i )f (x (k) j , z (k)
-j ) -f 2 0 , which further allows again for the estimation of all scaled indices, second-order total sensitivity indices and second-order Sobol indices. Although not mentioned there, it is obvious that the same applies for the total indices estimated by the Jansen formula instead of pick-freeze, which needs the same function evaluations.

Properties of Monte Carlo estimators

To assess the error coming from the Monte Carlo estimation, a straightforward idea is to use bootstrap [START_REF] Archer | Sensitivity measures, ANOVAlike techniques and the use of bootstrap[END_REF]. To do this, in the Monte Carlo estimation of an index of interest, the n summands in formula (2.14) are resampled with replacement for a high number of times (≈ 10 000). Each time, the corresponding index is computed and confidence intervals can be constructed from the resulting bootstrap sample. No new evaluations of f are necessary, as only the already present evaluations are used.

On the other hand, some theoretical results are sometimes available. Janon et al. n , uses the straightforward expressions f 0 (2.17) and D (2.18),

S 1 n = 1 n n k=1 f (x (k) )f (x (k) I , z (k) -I ) -1 n n k=1 f (x (k) ) 2 1 n n k=1 f (x (k) ) 2 -1 n n k=1 f (x (k) ) 2 , ( 2
S 2 n = 1 n n k=1 f x (k) f x (k) I , z (k) -I -1 n n k=1 f (x (k) )+f x (k) I ,z (k) -I 2 2 1 n n k=1 f (x (k) ) 2 +f x (k) I ,z (k) -I 2 2 -1 n n k=1 f (x (k) )+f x (k) I ,z (k) -I 2 2 .
(2.25) [START_REF] Janon | Asymptotic normality and efficiency of two sobol index estimators[END_REF] show that both estimators are consistent,

S 1 n a.s. -→ n→∞ D C I D , S 2 n a.s. -→ n→∞ D C I D ,
and that they are asymptotically normal in the form of

√ n S 1 n - D C I D d -→ n→∞ N 0, σ 2 1 , √ n S 2 n - D C I D d -→ n→∞ N 0, σ 2 2 .
See Janon et al. (2013, Prop. 3.3) for the specification of the variances σ 2 1 and σ 2 2 .

In their main result, they show that S 2 n is asymptotically efficient for the estimation of van der Vaart (1998, Chapter 25) in a framework with exchangeable variables.

D C I /D in the notion of

Frequency-based estimation

Beside direct Monte Carlo estimation, there is a further class of estimation methods, the frequency-based estimation, which is computationally cheaper than Monte Carlo estimation, but also slightly biased.

The first method introduced is the Fourier amplitude sensitivity test (FAST) by [START_REF] Cukier | Nonlinear sensitivity analysis of multiparameter model systems[END_REF], which allows for the estimation of first-order Sobol indices. Sample points of X are chosen such that the indices can be interpreted as amplitudes obtained by Fourier analysis of the function. More precisely, the design of n FAST points is such that

x (k) i = G i (sin(ω i s k )), i = 1, . . . , d, k = 1, . . . , n FAST , s k = 2π(k -1) n FAST ,
with G i functions to ensure that the sample points follow the distribution of X. The set of integer frequencies {ω 1 , . . . , ω d }, associated to the input variables, is chosen to be as free of interferences as possible, where free of interferences up to the order M means that p i=1 a i ω i = 0 for p i=1 |a i | ≤ M + 1 [START_REF] Tissot | Bias correction for the estimation of sensitivity indices based on random balance designs[END_REF]. In practice M = 4 or M = 6 are used.

For a weight ω, the Fourier coefficients for each input variable can then be numerically estimated by

A ω = 1 n FAST n FAST j=1 f (x(s j )) cos(ωs j ), B ω = 1 n FAST n FAST j=1 f (x(s j )) sin(ωs j ),
(2.26) and the first-order Sobol indices can be estimated by the sum of the corresponding amplitudes up to the order M ,

FAST D i = 2 M p=1 ( A 2 pω i + B 2 pω i ).
(2.27)

An estimate of the overall variance is given by the sum of all amplitudes,

FAST D = 2 n FAST /2 n=1 ( A 2 n + B 2 n ).
(2.28)

Extended FAST (EFAST), first presented in [START_REF] Saltelli | A quantitative model-independent method for global sensitivity analysis of model output[END_REF], is a method to also compute total sensitivity indices for single input variables using FAST. Here, a high frequency ω i is assigned to the variable in question X i and low frequencies, e.g. ω -i = 1, are assigned to all other variables. Then the total sensitivity index of X i is estimated over the amplitude for the low frequencies,

EFAST D T i = D -2 M p=1 ( A 2 pω -i + B 2 pω -i ).
A similar method, which avoids the problem of interferences, is RBD-FAST, where the RBD stands for random balanced design. RBD-FAST is a group of modifications of FAST, which use random permutations of design points to avoid interferences. The original idea by [START_REF] Tarantola | Random balance designs for the estimation of first order global sensitivity indices[END_REF] is to assign the same frequency to all input variables, but to randomize their values independently before evaluating f . The firstorder Sobol index of X i can then be estimated by reordering the evaluations in the way X i was permuted, so that the amplitude at the frequency returns the sensitivity of X i only. This works because the order of the evaluations is important in (2.26). Thus, all first-order Sobol indices can be estimated simultaneously with one set of evaluations.

In the same work a further RBD-FAST method called Hybrid is mentioned by which second-order and second-order closed sensitivity indices can be estimated. [START_REF] Mara | Extension of the RBD-FAST method to the computation of global sensitivity indices[END_REF] suggests a RBD-FAST method to compute the total sensitivity index of a group of variables D T I . Here, simple frequencies like ω = {1, . . . , d} are assigned to the variables. Then n RBD = 2(M d + L) design points are generated over a periodic curve, with L > 100 a tunable integer number, which regulates the sample size. The values of all variables in I are then randomly permuted, either differently per variable or identically, and the experiment is evaluated at the points. The total sensitivity index is estimated by

RBD D T I = n RBD L n RBD /2 p=dM +1 ( A 2 p + B 2 p ).
(2.29)

The EASI algorithm by Plischke ( 2010) is a further frequency-based method for the estimation of first-order Sobol indices. It has the advantage that no specific sampling scheme is necessary, which makes it possible to use already available data. It can be seen as a reverse RBD-FAST approach. The evaluations are reordered so that the x i are approximately triangular-shaped, which corresponds to the frequency ω = 1. Then the first-order Sobol index can be estimated as for RBD-FAST.

Derivative-based methods

This section shows different sensitivity analysis methods that, in contrast to variancebased methods, use the partial derivative of an input variable as an indicator of its influence. As the derivative quantifies the change in the output when the variable is changed, it can be used as a sensitivity measure.

Morris screening by elementary effects

In the field of input variable screening in computer experiments, that is finding the important variables out of a large number of variables using a minimal number of runs, Morris screening is one of the most popular methods [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF]. Among the advantages of the method over common screening methods are that it covers the entire input space and that it is not restricted to linear influences. Partial differences at different points in the input space, so-called elementary effects, are computed

d i (x) = |y(x 1 , . . . , x i + ∆ p , . . . , x k ) -y(x)| ∆ p ,
where the values of x are sampled on a grid with p levels for each input and ∆ p is a predetermined multiple of 1/(p -1). The distribution of the d i provides information on the behavior of X i . The mean, the so-called Morris Importance Measure, indicates the overall influence. The standard deviation indicates the linearity of the influence. A value close to 0 suggests linear behavior, a high value nonlinear or interaction behavior.

To obtain the distribution of elementary effects, the values of x are sampled randomly or by a more elaborate design that uses trajectories in the space of the input variables, starting from random points on the grid. This reduces the number of necessary runs almost by half compared to random sampling, resulting in a total number of r(d + 1) runs, with r the number of elementary effects per variable.

DGSM

A more accurate method based on Monte Carlo integration was introduced by [START_REF] Kucherenko | Monte Carlo evaluation of derivative-based global sensitivity measures[END_REF]. They use indices they call derivate-based global sensitivity measures (DGSM), first presented by [START_REF] Sobol | On an alternative global sensitivity estimator[END_REF] as an interpretation of variance-based indices. In contrast to Morris screening, the square of the partial derivative is used,

ν i = ∂f (X) ∂X i 2 dµ(X).
(2.30)

Like variance-based indices, the DGSM quantify the influence of the input variables and they share with the total sensitivity indices the ability that, (under the mild assumption that µ is continuous and its support is equal to the domain ∆ of X), if ν i = 0, then f (x) does not depend on x i . Indeed, [START_REF] Sobol | Derivative based global sensitivity measures and their link with global sensitivity indices[END_REF] could show a connection between the two sensitivity methods for the uniform and normal distributions, which was extended by [START_REF] Lamboni | Derivative-based global sensitivity measures: general links with Sobol' indices and numerical tests[END_REF] for general continuous distributions. It states that (2.31) provided that µ belongs to a class of distributions that satisfy a Poincaré inequality,

D T i ≤ C(µ i )ν i ,
g(X) 2 dµ(X) ≤ C(µ) ∇g(X) 2 dµ(X), (2.32)
for all functions g in L 2 (µ) such that g(X) dµ(X) = 0, and ∇g ∈ L 2 (µ). A

Poincaré constant C(.) is a characteristic of the corresponding distribution µ. The best possible constant for a given distribution, the optimal Poincaré constant, is denoted by C opt (µ). When C(µ) = C opt (µ), there can exist certain functions f opt for which the Poincaré inequality is an equality. See [START_REF] Lamboni | Derivative-based global sensitivity measures: general links with Sobol' indices and numerical tests[END_REF] or [START_REF] Roustant | Crossed-derivative based sensitivity measures for interaction screening[END_REF] for a comprehensive summary.

Formula (2.31) implies that the DGSM can be used as upper boundaries for the total sensitivity indices and thus can serve as a cheaper method for the screening of noninfluential input variables since, after [START_REF] Lamboni | Derivative-based global sensitivity measures: general links with Sobol' indices and numerical tests[END_REF], the DGSM seems to be computationally more tractable than variance-based measures, especially for problems with higher numbers of inputs. This applies especially when the computer experiment additionally provides the partial derivatives. Then, the Monte Carlo integration can be performed directly. If not, the derivatives can be approximated numerically by finite differences. The DGSM for a n × d data set x and a small real number δ * can be estimated by

ν i = 1 n n k=1 f (x (k) i + δ * , x (k) -{i} ) -f (x (k) ) δ * 2 .
To explore the difference between the DGSM and variance-based indices, let us compare two one-dimensional case functions,

f 1 (x) = 1 1 + 49 exp(-2 log(49)x) , f 2 (x) = 1 2 sin(50x) + 1 2 .
(2.33)

Plots are shown in Fig. 2.1. In this one-dimensional situation, the total index, the main index, and the overall variance coincide, the variance-based index reduces to

D 1 = D T 1 = D = 1 0 f (x) 2 dx -   1 0 f (x) dx   2 ,
and for the DGSM we have

ν 1 = 1 0 (f (x)) 2 dx.
When computing these values for the two functions, we see that the DGSM is much higher for function f 2 (≈ 310.9) than for f 1 (≈ 1.3), as it visibly varies much more.

Nevertheless, the overall variance is for both functions almost the same (≈ 0.126), as the function evaluations over the space of x share a similar variance. The variancebased index ignores the local variance. This highlights the difference of the notion sensitivity for variance-based indices and DGSM. In an hypothetical underlying model y = f 1 (x 1 )+f 2 (x 2 ), the DGSM would rank the second variable much higher whereas the 0.0 0.2 0.4 0.6 0.8 1.0 0.0 0.2 0.4 0.6 0.8 1.0

x f 1 (x), f 2 (x) f 1 (x) f 2 (x)

Ongoing developments in sensitivity analysis

To complete the overview of methods in sensitivity analysis, we want to show some interesting topics of current research in order to give a little insight into the progression of the field.

Owen (2013b) creates a general framework for variance-based indices by using linear combinations of Sobol indices. The various possible estimators for the different indices are summarized and categorized, which leads to some new estimators. Among others, it supplies a bias-corrected pick-freeze estimator and estimators for Sobol indices with a reduced number of necessary function evaluations.

For the case of multivariate output, several ideas exist that incorporate correlations between output variables to global sensitivity analysis. Indices were introduced by [START_REF] Lamboni | Multivariate sensitivity analysis to measure global contribution of input factors in dynamic models[END_REF] and studied further by [START_REF] Gamboa | Sensitivity indices for multivariate outputs[END_REF] that extend sensitivity indices to multivariate output, which -relating to the FANOVA decomposition -base on the decomposition of the covariance of the output. An overview of different approaches to the analysis of multivariate output can be found in [START_REF] Garcia-Cabrejo | Global sensitivity analysis for multivariate output using polynomial chaos expansion[END_REF].

Control variates is a technique to reduce the variance in Monte Carlo integration.

Applying control variates to the estimation of total sensitivity indices via the Jansen formula, [START_REF] Kucherenko | Application of the control variate technique to estimation of total sensitivity indices[END_REF] find the formula

D T i = 1 2 E [f (X) -f i (X i ) -(f (Z i , X -i ) -f i (Z i ))] 2 + D i ,
which can indeed improve the estimators efficiency. It requires the FANOVA decomposition term f i as well as the first-order Sobol index D i . If not known analytically, these terms have to be extracted from metamodels.

A further interesting idea is an extension to Sobol indices that takes the goal of further analysis into consideration [START_REF] Fort | New sensitivity analysis subordinated to a contrast[END_REF]. The idea is that the estimation of a mean or a median could involve different input variables than the estimation of extreme quantiles. They set up a framework they call Goal Oriented Sensitivity Analysis, where new sensitivity indices for each statistical purpose are defined by applying contrast functions. When considering the estimation of α-quantiles, a function of indices over α is returned.

This chapter presented the basic sensitivity analysis techniques used for the exploration of input variables in computer experiments. The following chapter will present and examine methods that focus on the exploration of interactions between variables.

Chapter 3

Sensitivity Analysis for Interaction Screening

The total interaction index (TII) is a sensitivity index for interaction screening in computer experiments, which can further be used for visualization and block-additive decomposition of the computer model. As the total sensitivity index, introduced in Chapter 2, is ideal for input screening, since it sums up the total influence of the input variables, the screening of interactions can be done by the TII. It contains, for any pair of input variables, the influence of the interaction of the pair plus all interactions containing the pair. Thus, corresponding to the total sensitivity index where a value close to zero is a strong indicator that a variable can be excluded [START_REF] Saltelli | Sensitivity analysis practices: Strategies for model-based inference[END_REF], a TII close to zero indicates that the two variables do not interact. It has been used but not investigated closer in [START_REF] Mühlenstädt | Data-driven Kriging models based on FANOVA-decomposition[END_REF] in order to set up so-called FANOVA graphs and Kriging models with block-additive kernels. [START_REF] Liu | Estimating mean dimensionality of analysis of variance decompositions[END_REF] introduced a more general index for uniform distributions as a measure of importance of interactions, which was also used in [START_REF] Hooker | Discovering additive structure in black box functions[END_REF] in the data mining framework.

This chapter starts with an example that gives a direct motivation to the TII. Then, several variance-based estimation methods for the TII are introduced. Their properties are analyzed theoretically as well as on simulations. The problem of finding a threshold is addressed as well as the implementation within the statistical software programming 25 environment R (R Core Team, 2014). In addition, crossed DGSM, extensions of the DGSM to second-order interactions, are presented, which provide upper bounds for the TII.

Some results from this chapter are published in the contributions "Total interaction index: a variance-based sensitivity index for second-order interaction screening" (Fruth et al., 2014b) and "Crossed-derivative based sensitivity measures for interaction screening" [START_REF] Roustant | Crossed-derivative based sensitivity measures for interaction screening[END_REF].

Motivational example

With the TII, it is possible to discover the block-additive structure of the underlying function, that is we can identify a decomposition into groups of input variables such that variables between groups do not interact. As an illustration, we analyze the following function,

f (X 1 , . . . , X 6 ) = cos([1, X 1 , X 3 , X 5 ] α) + sin([1, X 2 , X 4 , X 6 ] γ), with X i i.i.d. ∼ U [-1, 1], i = 1, . . . , 6, α = [-0.8, -1.1, 1, 1.1] and γ = [-0.51, 0.9, -1.1] ,
where the prime stands for the transpose.

Clearly, the variables in the group {X 1 , X 3 , X 5 } do not interact with the variables in the group {X 2 , X 4 , X 6 }, which induces a block-additive structure of the form f (x) =

f 135 (x 1 , x 3 , x 5 ) + f 2,4,6 (x 2 , x 4 , x 6
). This, however, cannot be seen from the common first-order and total sensitivity indices (Fig. 3.1, left). We only see that all variables have first-order and total effects, but do not know how variables interact with each other.

If we now estimate the TII for each combination of input variables, we obtain such information, which can be conveniently plotted in a so-called FANOVA graph in Fig. 3.1 on the right [START_REF] Mühlenstädt | Data-driven Kriging models based on FANOVA-decomposition[END_REF]. The TII of the two variables is represented by the thickness of the edge between the two corresponding vertices. With the graph, the partition into the two additive groups is clearly visible.

X 1 X 2 X 3 X 4 X 5 X
This information about the block-additive interaction structure of a function can be exploited in several fields of computer experiment analysis. [START_REF] Mühlenstädt | Data-driven Kriging models based on FANOVA-decomposition[END_REF] show that it can be used to improve Kriging model predictions by adapting the Kriging kernel to the block-additive structure. For the example above that would mean to modify the kernel k from Equation (2.3) as follows,

k(h) = k 1 (h 2 , h 4 , h 6 ) + k 2 (h 1 , h 3 , h 5 ).
Furthermore, the block-additive structure can be exploited to simplify and parallelize optimization, as described in [START_REF] Ivanov | A parallel optimization algorithm based on FANOVA decomposition[END_REF]. There, the idea is to optimize the separate groups independently, which reduces the optimization dimensions and enables parallelization, an important property as optimization techniques are usually sequential. For the example above, the six-dimensional optimization problem of minimizing f simplifies into two three-dimensional ones, where the c 1 , . . . , c 6 are constants, which are not varied in the optimization, min

x 1 ,...,x 6 f (x 1 , . . . x 6 ) = min x 1 ,x 3 ,x 5 f (x 1 , c 2 , x 3 , c 4 , x 5 , c 6 ) + min x 2 ,x 4 ,x 6 f (c 1 , x 2 , c 3 , x 4 , c 5 , x 6 ).

Total interaction indices

The TII measures the total influence of a second-order interaction between two input variables. It contains, for any pair of input variables, the influence of their interaction plus all interactions containing both indices. This is different from the second-order Sobol index (2.8), which does not contain higher interactions, and from the secondorder total sensitivity index (2.10), which additionally contains first-order effects as well as higher-order interactions of only one of the two variables.

Definition 3.1. The total interaction index (TII) of two input variables X i and X j is defined by

D i,j = Var   I ⊇{i,j} f I (X I )   = I ⊇{i,j} D I . (3.1)
The TII is a special second-order version of the more general superset importance, which was introduced by [START_REF] Liu | Estimating mean dimensionality of analysis of variance decompositions[END_REF] for uniform distributions (Υ 2 u ) as a measure of importance of interactions and their supersets. It was also investigated by [START_REF] Hooker | Discovering additive structure in black box functions[END_REF] in the data mining framework (σ 2 u ). The superset importance Υ 2 u is defined for any subset u ⊆ {1, . . . , d} as

Υ 2 u = σ 2 u = I⊇u D I . (3.2)
One important property of the TII is its link to total sensitivity indices and closed indices respectively.

Proposition 3.1.

D i,j = D T i + D T j -D T i,j , (3.3)
D i,j = D + D C -{i,j} -D C -i -D C -j . (3.4)
Proof. The proof of (3.3) is obtained by simple set transformation. Since

I⊇{i}∨I⊇{j} D I = I⊇{i} D I + I⊇{j} D I - I⊇{i,j} D I , it holds that I⊇{i,j} D I = I⊇{i} D I + I⊇{j} D I - I⊇{i}∨I⊇{j} D I .
Equation (3.4) then can be deduced from (3.3) using the connection between total and closed indices (2.13),

D i,j = D T i + D T j -D T i,j = (D -D C -i ) + (D -D C -j ) -(D -D C -{i,j} ) = D + D C -{i,j} -D C -i -D C -j .
From the connection to the superset importance, another representation of the TII is given by [START_REF] Liu | Estimating mean dimensionality of analysis of variance decompositions[END_REF].

Proposition 3.2. Liu and Owen's formula

D i,j = 1 4 E f (X i , X j , X -{i,j} ) -f (X i , Z j , X -{i,j} ) -f (Z i , X j , X -{i,j} ) + f (Z i , Z j , X -{i,j} ) 2 , (3.5)
where Z i (resp. Z j ) is an independent copy of X i (resp. X j ).

Proof. A proof is given in [START_REF] Liu | Estimating mean dimensionality of analysis of variance decompositions[END_REF]. In addition, (3.5) can be connected to (3.4). When expanding the squared sum in (3.5), of the 10 resulting terms the four squared terms are simply equal to E(f (X) 2 ) = D + f 2 0 . The six double products gather two by two, and result in the pick-freeze formula (2.15)

• E[f (Z i , X j , X -{i,j} )f (Z i , Z j , X -{i,j} )] = E[f (X i , X j , X -{i,j} )f (X i , Z j , X -{i,j} )] = D C -j + f 2 0 • E[f (X i , Z j , X -{i,j} )f (Z i , Z j , X -{i,j} )] = E[f (X i , X j , X -{i,j} )f (Z i , X j , X -{i,j} )] = D C -i + f 2 0 • E[f (Z i , X j , X -{i,j} )f (X i , Z j , X -{i,j} )] = E[f (X i , X j , X -{i,j} )f (Z i , Z j , X -{i,j} )] = D C -{i,j} + f 2 0
Finally, combining all the terms gives (3.4).

Remark: Connection to ANOVA [START_REF] Liu | Estimating mean dimensionality of analysis of variance decompositions[END_REF] state that formula (3.5) "can also be obtained through the classical formulas for expected mean squares in the discrete ANOVA as established by [START_REF] Cornfield | Average values of mean squares in factorials[END_REF]". This connection becomes evident when looking at the formula for the interaction Sum of Squares in an ANOVA with two input variables A and B (Cornfield and Tukey, 1956, p. 916),

SS AB = n a i=1 b j=1 (ȳ ij• -ȳi•• -ȳ•j• + ȳ••• ) 2 ,
with a and b the number of levels and y ijk , i = 1, . . . , a, j = 1, . . . , b, k = 1, . . . , n the realized output values with n the number of repetitions.

For a further TII computation notion, say we keep all variables other than X i and X j fixed and look at the Liu and Owen formula for the resulting two-dimensional function

f fixed (X i , X j ), D i,j = 1 4 E (f fixed (X i , X j ) -f fixed (X i , Z j ) -f fixed (Z i , X j ) + f f ixed (Z i , Z j )) 2 .
Then we see that the TII of the original function can be computed as the expected value over this second-order interaction index of the two-dimensional function obtained by fixing all variables except X i and X j . Thus, the Liu and Owen formula is a so-called fixing method, introduced by [START_REF] Mühlenstädt | Data-driven Kriging models based on FANOVA-decomposition[END_REF].

Proposition 3.3. (Fixing method). For any x -{i,j} , define f fixed as the two-dimensional function f fixed : (x i , x j ) → f (x) obtained from f by fixing all input variables except x i and x j . Let D i,j|X -{i,j} denote the second-order Sobol index of f fixed (X i , X j ), which depends on the fixed variables X -{i,j} . Then the TII of X i and X j is obtained by integrating D i,j|X -{i,j} with respect to X -{i,j} ,

D i,j = E D i,j|X -{i,j} . (3.6)
Proof. Since the function f fixed is two-dimensional, it has only one interaction, which is a second-order one, and coincides with its TII. Hence, this interaction can be computed by applying (3.5) to f fixed ,

D i,j|X -{i,j} = 1 4 E [f fixed (X i , X j ) -f fixed (X i , Z j ) -f fixed (Z i , X j ) + f fixed (Z i , Z j )] 2 .
Now we can rewrite the right hand side by using conditional expectations,

D i,j|X -{i,j} = 1 4 E f (X i , X j , X -{i,j} ) -f (X i , Z j , X -{i,j} ) -f (Z i , X j , X -{i,j} ) + f (Z i , Z j , X -{i,j} ) 2 |X -{i,j} .
Taking the expectation with respect to X -{i,j} gives the result.

TII estimation

The different representations of the TII derived in the last section and the estimation methods of the common indices from Section 2.3 are now combined to construct estimation methods for the TII, followed by remarks on their properties.

Estimation via total sensitivity indices

To apply the link between the TII and total sensitivity indices from Prop. 3.1, formula (3.3), estimators for the total sensitivity index of first-and second-order are required.

Thus, either the Monte Carlo estimator via the Jansen formula (2.23) or the RBD-FAST method (2.29) can be used, resulting in the TII estimators Jansen TII estimator and RBD-FAST TII estimator,

Jan D i,j = Jan D T i + Jan D T j - Jan D T {i,j} (3.7) RBD D i,j = RBD D T i + RBD D T j - RBD D T {i,j} (3.8)
For the Jansen TII estimator, evaluations can be saved by computing the last term Jan D T {i,j} via f (x -i , z i ) and f (x -j , z j ), which have been evaluated in the computation of the previous terms, rather than via f (x) and f (x -{i,j} , z {i,j} ). The approach is discussed further in Section 3.4.

Estimation via closed sensitivity indices

The second link from Prop. 3.1 between the TII and the closed sensitivity indices (3.4) requires the estimation of high-order closed indices D C -i and D C i,j . To this point, no reliable FAST or RBD-FAST method exists, so Monte Carlo integration is considered.

As the required closed indices are of high order, they are expected to be large. Thus, the estimator based on the pick-freeze formula (2.15, 2.16) is more suitable here than the one based on the strategies from Correlation 1 and Correlation 2 (2.21) since they are not recommended when the index to be estimated is rather large (Owen, 2013a).

The pick-freeze TII estimator of the TII is constructed as

pf D i,j = D + pf D C -{i,j} - pf D C -i - pf D C -j .
(3.9)

As for the Jansen TII estimator, evaluations can be reused in the second-order term by using f (x -i , z i ) and f (x -j , z j ) rather than f (x) and f (x -{i,j} , z {i,j} ).

Estimation via Liu and Owen's formula

Another TII estimation method is directly obtained by Liu and Owen's formula in Prop. 3.2, the Liu and Owen TII estimator

LO D i,j = 1 4 × 1 n n k=1 f (x k i , x k j , x k -{i,j} ) -f (x k i , z k j , x k -{i,j} ) -f (z k i , x k j , x k -{i,j} ) + f (z k i , z k j , x k -{i,j} ) 2 . (3.10)
Estimation via fixing method Following Prop. 3.3, the TII can be computed by averaging the second-order interaction of two-dimensional functions in the way it has been done in [START_REF] Mühlenstädt | Data-driven Kriging models based on FANOVA-decomposition[END_REF] by the following scheme:

Let us consider a couple of integers (i, j), with i < j. For k = 1, . . . , n MC carry out the following steps.

1. Simulate x k -{i,j} from the distribution of X -{i,j} , that is take a single sample of all input variables except X i and X j .

Create the two-dimensional function

f f ixed by fixing f on x k -{i,j} f fixed (X i , X j ) = f (x k 1 , . . . , X i , . . . , X j , . . . , x k d ).
3. Using the FAST estimator (Section 2.3.3), compute the second-order Sobol index

of f fixed , denoted D k i,j|X -{i,j}
, by

D k i,j|X -{i,j} = FAST D k |X -{i,j} - FAST D k i|X -{i,j} - FAST D k j|X -{i,j} , TII estimator derived from formula Jansen Prop. 3.1 (3.3) Jan D i,j = Jan D T i + Jan D T j - Jan D T {i,j} RBD-FAST Prop. 3.1 (3.3) RBD D i,j = RBD D T i + RBD D T j - RBD D T {i,j} pick-freeze Prop. 3.1 (3.4) pf D i,j = pf D + pf D C -{i,j} - pf D C -i - pf D C -j
Liu and Owen Prop. 3.2 the first-order indices of X i and X j , respectively.

LO D i,j = 1 4 × 1 n n k=1 f (x k i , x k j , x k -{i,j} ) -f (x k i , z k j , x k -{i,j} ) -f (z k i , x k j , x k -{i,j} ) + f (z k i , z k j , x k -{i,j} ) 2 fixing method Prop. 3.3 fix D i,j = 1 n MC n MC k=1 D k i,j|X -{i,j}
After carrying out the three points for all k = 1, . . . , n MC , compute the fixing method TII estimator by

fix D i,j = 1 n MC n MC k=1 D k i,j|X -{i,j}
.

(3.11)

A summary of all considered TII estimators can be found in Tab. 3.1.

Theoretical properties of TII estimators

In the following, the introduced TII estimators are compared according to different statistical properties.

Nonnegativity

As the indices measure the sum of variances, the estimates should be nonnegative.

This holds for the Liu and Owen TII estimator (3.10) which is a sum of squares.

However, negative estimates can occur for the Jansen (3.7), the pick-freeze (3.9) and the RBD-FAST TII estimator (3.8). For the fixing method TII estimator (3.11), there is a sufficient condition, which results from the following proposition:

Proposition 3.4. Let f be a two-dimensional function, and consider its second-order

interaction D 12 = D -D 1 -D 2 . Denote by D 12 = FAST D - FAST D 1 - FAST D 2 its FAST
estimate (see Section 2.3.3) and assume that:

(i) ω 1 and ω 2 are free of interference up to order 2M ,

(ii) N ≥ 2M × max(ω 1 , ω 2 ).
Then D 12 ≥ 0.

Proof. Denote the sets W ω i ,M = {pω i , p = 1, . . . , M } for i = 1, 2, and W N = {1, . . . , N/2}. With (2.27) and (2.28), we have

D 12 /2 = n∈W N ( A 2 n + B 2 n ) - n∈W ω 1 ,M ( A 2 n + B 2 n ) - n∈W ω 2 ,M ( A 2 n + B 2 n ). Now, the condition (i) ensures that W ω 1 ,M ∩ W ω 2 ,M = ∅ while (ii) implies that W ω i ,M ⊆ W N , for i = 1, 2. Hence, D 12 /2 = n∈W N -(W ω 1 ,M ∪W ω 2 ,M ) ( A 2 n + B 2 n ) ≥ 0.
It is a direct consequence of Prop. 3.4 that if (i) and (ii) are satisfied, then (3.11) returns positive values. In practice, one can use for instance ω 1 = 11, ω 2 = 35 [START_REF] Mara | Extension of the RBD-FAST method to the computation of global sensitivity indices[END_REF], which are free of interferences up to 2M for the usual orders M = 4, 6. Then the minimal value of N is 2 × 6 × max{11, 35} = 420.

Bias

The methods based on direct Monte Carlo estimation, Jansen and Liu and Owen TII estimator, are unbiased since only direct mean estimators are used for the conditional expectations. This is especially remarkable in combination with the positivity of the Liu and Owen TII estimator, as it implies that when the true value is zero, the estimator is identical to zero as well.

Proposition 3.5. If D i,j = 0, then the Liu and Owen TII estimator is equal to zero:

LO D i,j ≡ 0.
Proof. Starting with the FANOVA decomposition (2.4), if D i,j = 0, then all of the terms containing both x i and x j vanish. So the decomposition reduces to

f (x) = i / ∈I∧j / ∈I f I (x I ) + i∈I∧j / ∈I f I (x I ) + i / ∈I∧j∈I f I (x I ) = a 0 (x -{i,j} ) + a i (x i , x -{i,j} ) + a j (x j , x -{i,j} ).
Inserting this representation of f into the squared term of the Liu and Owen TII estimator (3.10) returns zero,

f (x k i , x k j , x k -{i,j} ) -f (x k i , z k j , x k -{i,j} ) -f (z k i , x k j , x k -{i,j} ) + f (z k i , z k j , x k -{i,j} ) = a i (x k i , x k -{i,j} ) -a i (x k i , x k -{i,j} ) -a i (z k i , x k -{i,j} ) + a i (z k i , x k -{i,j} ) + a j (x k j , x k -{i,j} ) -a j (z k j , x k -{i,j} ) -a j (x k j , x k -{i,j} ) + a j (z k j , x k -{i,j} ) = 0.
In the pick-freeze TII estimator, the estimation of f 2 0 is necessary in the estimation of the closed sensitivity indices (2.16). If we estimate f 2 0 by taking the square of f 0 (2.17) or * f 0 (2.19), we introduce a bias of -Var(f 0 ), as was noted by Owen (2013b),

[E(f 0 )] 2 = E (f 0 ) 2 -Var(f 0 ).
However, as long as E(f (X) 4 ) < ∞, the bias is asymptotically negligible. Owen (2013b) remarks that the bias may be important in the case of small closed sensitivity indices, but those are not expected in the estimation of pf D i,j , as mentioned in Section 3.2. To nonetheless ensure unbiasedness, unbiased estimators for the closed sensitivity index can be used such as Frequency-based estimators are generally prone to bias [START_REF] Tissot | Bias correction for the estimation of sensitivity indices based on random balance designs[END_REF]. This bias might even be enhanced here through the use of a combination of FAST estimators for the fixing method TII estimator and RBD-FAST estimators for the RBD-FAST TII estimator.

Asymptotic properties of the Liu and Owen TII estimator

Corresponding to the asymptotic properties of the pick-freeze estimator of closed sensitivity indices (end of Section 2.3.2), it is possible to derive asymptotic properties of the Liu and Owen TII estimator (3.10). To do this, the estimator for a pair of input variables {X i , X j } is written as

T n = LO D i,j = 1 n n k=1 ∆ k i,j 2 4 , with ∆ k i,j = f (X k i , X k j , X k -{i,j} ) -f (X k i , Z k j , X k -{i,j} ) -f (Z k i , X k j , X k -{i,j} ) + f (Z k i , Z k j , X k -{i,j} ).
Proposition 3.6. Let P be the set of all cumulative distribution functions of exchangeable random vectors in L 2 (R 2 ), that is for a P ∈ P it holds for any random vectors X and X that P (X, X ) = P (X , X). Then the following propositions hold for the Liu and Owen TII estimator

T n . a) T n is consistent for D i,j , T n a.s. -→ n→∞ D i,j . b) T n is asymptotically normally distributed, √ n (T n -D i,j ) d -→ n→∞ N 0, Var[(∆ 1 i,j ) 2 ] 16 .
c) T n is asymptotically efficient in the notion of van der Vaart (1998) for estimating D i,j for P ∈ P.

Proof. The results a) and b) are a direct application of the law of large numbers and the central limit theorem, applied to the variables (∆ k i,j ) 2 .

Result c) follows from the fact that estimators with symmetrical expressions are asymptotically efficient in the framework of exchangeable variables, as proved by Janon et al.

(2013) in their Lemma 2.6 (2):

Let Φ 2 : R 2 → R be a symmetric function in L 2 (P ) and f a deterministic function on R ⊂ R p 1 +p 2 of independent random input variables X ∈ R p 1 and Z ∈ R p 2 . Z denotes an independent copy of Z and X k , Z k , Z k , k = 1, . . . , n denote independent samples of the corresponding variables. The sequence {Φ 2 n } n∈N given by

Φ 2 n = 1 n n k=1 Φ 2 f (X k , Z k ), f (X k , Z k ) is asymptotically efficient for estimating E(Φ 2 (f (X, Z), f (X, Z ))) for P ∈ P. More precisely, denote X k = (X k j , Z k j , X k -{i,j} ), Z k = X k i , Z k = Z k i , and let g be the function defined over R d × R by g(a, b) = f (b, a 1 , a 3 , . . . , a d ) -f (b, a 2 , a 3 , . . . , a d ). Then ∆ k i,j = g(X k , Z k ) -g(X k , Z k ).
Therefore

T n = 1 n n k=1 Φ 2 (g(X k , Z k ), g(X k , Z k )),
and

D i,j = E(Φ 2 (g(X 1 , Z 1 ), g(X 1 , Z 1 ))),
where Φ 2 is the two-dimensional function defined over R 2

Φ 2 (u, v) = (u -v) 2 4 .
Remark that Z k and Z k are independent copies of each other, both independent of X k , and that Φ 2 is a symmetric function. With the following change in notation

i ← k, X ← X , Z ← Z, Z ← Z , f ← g
the result follows from the Lemma of Janon et al. ( 2013) named above.

The two last propositions can be extended to the general superset importance (3.2), including the case of the total sensitivity index of one input variable.

Proposition 3.7. Let Υ I = J⊇I D J be the superset importance for a set I. Define Then T I,n is asymptotically normal and asymptotically efficient for Υ I .

T I,n = 1 n n k=1 (∆ k I ) 2 2 |I| , with ∆ k I = J⊆I (-1) |I-J| f (Z k J , X k -J ),
Proof. Note that T I,n is the sample version of the formula (10) given by [START_REF] Liu | Estimating mean dimensionality of analysis of variance decompositions[END_REF] for Υ I (with suitable change of notations). The proof of asymptotic normality is thus a direct consequence of the central limit theorem. For asymptotic efficiency, the proof relies on arguments similar to Prop. 3.6.

• When I = {i} is a single input variable, we have

∆ k I = f (Z k i , X k -i ) -f (X k i , X k -i ), which is of the form g(X k , Z k ) -g(X k , Z k ) with Z k = Z k i , Z k = X k i , X k = X k -i , and g(a, b) = f (b, a).
• When |I| ≥ 2, let us choose i ∈ I. Then, by splitting the subsets of I into two parts, depending whether they contain {i}, we have

∆ k I = J⊆I-{i} (-1) |I-J| f (Z k J , X k -J ) + J⊆I-{i} (-1) |I-(J∪{i})| f (Z k J∪{i} , X k -(J∪{i}) ) = J⊆I-{i} (-1) |I-J| f (Z k J , X k i , X k -(J∪{i}) ) - J⊆I-{i} (-1) |I-J| f (Z k J , Z k i , X k -(J∪{i}) ), which is also of the form g(X k , Z k ) -g(X k , Z k ) with Z k = X k i , Z k = Z k i , X k = (X k I-{i} , Z k I-{i} , X k -I
), and a suitable g since the second term in the difference is obtained from the first one by exchanging X k i and Z k i .

The result is then obtained by applying Lemma 2.6 in [START_REF] Janon | Asymptotic normality and efficiency of two sobol index estimators[END_REF] to the

symmetric function Φ 2 (u, v) = (u-v) 2
2 |I| , remarking that Z k and Z k are independent copies of each other, both independent of X k .

Corollary 3.1. It follows from Prop. 3.7 that the Jansen estimator of the total sensitivity index of a single variable D T i is asymptotically efficient.

Estimating the full set of TIIs

In most of the practical applications, the full set of d(d-1) 2 TIIs, that is the TII of each possible pair of the d input variables, is required simultaneously. The function evaluations, especially in the case without a metamodel, where f represents the actual computer experiment, are usually the most time-consuming part of the estimation.

Table 3.2 therefore gives an overview of the number of evaluations N that are required by each method for the estimation of the full set of TIIs, depending on the different parameter settings for each method.

In each of the n Monte Carlo runs, the Liu and Owen TII estimator requires one estimate corresponding to f (x), d estimates corresponding to f (z i , x -i ) -and at the same time to f (z j , x -j ) -and d 2 estimates corresponding to f (z i , z j , x -{i,j} ). Similarly, the RBD-FAST TII estimator requires d estimates of For the Jansen and the pick-freeze TII estimator, strategies that reuse computations corresponding to the strategy by [START_REF] Saltelli | Making best use of model evaluations to compute sensitivity indices[END_REF] (Section 2.3.2) can be applied, as addressed in their definitions (3.7) and (3.9). In the estimation of the first-and secondorder total indices for the Jansen TII estimator as well as in the estimation of the (d -1)-and (d -2)-order closed indices for the pick-freeze TII estimator, only the

(d + 1)n samples f (x k ) and the f (z k i , x k -i ), 1 ≤ k ≤ n, 1 ≤ i ≤ d are required
for the estimation of all TII indices. In addition, those function evaluations can be used to estimate the full set of first-order total sensitivity indices using either (2.22) or (2.23).

By adding only the n samples f (z k ), 1 ≤ k ≤ n, all first-order Sobol indices (2.16) can be estimated as well.

Table 3.2 clearly shows that in contrast to the other estimators, the Jansen as well as the pick-freeze TII estimator are scalable in the sense that the number of required runs increases linearly with the number of input variables d.

TII estimator number of function evaluations

Jansen N = (d + 1) × n RBD-FAST N = d 2 + d × 2(M d + L) pick-freeze N = (d + 1) × n Liu and Owen N = d 2 + d + 1 × n fixing method N = d 2 × n MC × n FAST

Simulation study

In the following, the behavior of the different estimators is investigated in simulations on known analytical real-valued functions. In order to cover a wide range of functions, functions are chosen that fall in different categories according to different properties based on [START_REF] Kucherenko | Monte Carlo evaluation of derivative-based global sensitivity measures[END_REF]. The complexity of the function is considered in two different categories, the interaction order and the functional form. For the latter, simple forms like low-order polynomial functions are considered as of low and complex nonlinear forms as of high complexity. In detail, the categories are

• the number of input variables: low (d ≤ 5) or high (d ≥ 10),

• interaction behavior: low (only low-order terms are dominant) or high (important high-order interactions),

• functional form: simple (simple functional form like low-order polynomial) or complex (complex nonlinear behavior). (1990), and Branin4 from [START_REF] Lehman | Designing computer experiments to determinate robust control variables[END_REF]. Table 3.4 shows how the functions are classified into the defined categories.

name function µ(X) function 1 sin(X 1 + X 2 ) + 0.4 cos(X 3 + X 4 ) U [-1, 1] 4 Ishigami sin(X 1 ) + 7 sin 2 (X 2 ) + 0.1X 4 3 sin(X 1 ) U [-π, π] 3 pure 3 rd order √ 3 3 3 i=1 X i U [-1, 1] 3 Branin4 1 30 f Branin (X 1 , X 2 )f Branin (X 3 , X 4 ) + (X 1 -π) 2 with f Branin (X 1 , X 2 ) = X 2 -5.1 4π 2 X 2 1 + 5 π X 1 -6 2 +10 1 -1 8π cos(X 1 ) + 10 X 1 , X 3 ∼ U [-5, 10], X 2 , X 4 ∼ U [0, 15]
high 2 nd order 3( For each of the functions, the TIIs are estimated a hundred times by all 5 methods. For each estimation, 2 000 Monte Carlo draws per index are used, 2 000× d 2 draws in total. The parameters of each method are determined using the connections in Tab. One reason for this might be the bias for RBD-FAST methods, mentioned in Section 3.3. This problematic behavior of the RBD-FAST TII estimator continues in the other four simulations. The other four estimators perform well in terms of bias, confirming the results of Section 3.3. The described bias of the pick-freeze TII estimator is not noticable, the deviations from the true mean are comparable to the Jansen and the Liu-Owen TII estimator.

7 i=1 X i )( 14 j=8 X j ) U [-1, 1] 14 high 14 th order √ 3 14 14 i=1 X i U [-1, 1] 14
Simulation boxplots for functions with high order interactions but a rather small number of input variables are shown in Fig. 3.3. Here, the Liu and Owen and the fixing method TII estimator show a stronger variance than before. The reason for this lies in the variance of the estimates D k i,j|X -{i,j} in the basic fixing method formula (3.6), which applies for both methods. For second-order interactions, those estimates do not depend on the fixed values X -{i,j} and thus vary only slightly while for higher interactions the estimates should differ with the fixed variables included in the interaction.

The much larger variance of the fixing method TII estimator can be additionally explained by the fact that the two steps, second-order interaction estimation and overall expectation, are performed separately, and the values for X -{i,j} are kept constant in each second-order index estimation, whereas in the Liu and Owen method both steps are performed simultaneously and thus are more efficient. Generally, for complex functions, the methods Jansen, Liu and Owen and pick-freeze TII estimator seem to perform equally well.

In the simulation study of the two high-dimensional test functions, a 10 times higher number, 20 000, of Monte Carlo draws is necessary due to numerical problems for the extremely complex case of a 14 th order interaction. As the RBD-FAST TII estimator This is due to the strategy to reuse computations, which makes those two methods scalable (see Section 3.4). The fixing method TII estimator shows a bias, perhaps caused by the bias in FAST described e.g. in [START_REF] Tissot | Bias correction for the estimation of sensitivity indices based on random balance designs[END_REF].

To get an impression of the asymptotic behavior, Fig. 3.5 shows exemplarily the estimates of the four reliable estimators Jansen, Liu and Owen, pick-freeze and fixing method TII estimator for index D 13 of the pure 3 rd order function. For each method, the index is estimated starting with the lowest number of evaluations possible and then evaluations are added step by step. The estimate of the asymptotically efficient estimator by Liu and Owen is the fastest to come close the true value followed by the pick-freeze and the Jansen TII estimator while the estimate of the fixing method TII estimator does not reach it within the 50 000 evaluations.

The simulation study shows that the Liu and Owen TII estimator indeed provides good results for almost all situations and can thus be recommended as method of choice. The only exception are situations with a high number of variables and strong interactions.

There, the pick-freeze TII estimator is rather preferable.

Threshold decision

Errors in the initial metamodel and in the index estimation usually perturb the TII estimation so that, when applying them to interaction screening, inactive indices are not necessarily zero. Therefore, a threshold needs to be introduced to separate active from inactive interactions. The decision rule is then to include the interaction if

D i,j D > δ.
(3.12)

If the TII estimation is further used for adapted modeling [START_REF] Mühlenstädt | Data-driven Kriging models based on FANOVA-decomposition[END_REF] or parallelized optimization [START_REF] Ivanov | A parallel optimization algorithm based on FANOVA decomposition[END_REF], it is more problematic to cut off an active interaction than to falsely accept an inactive interaction as active. For modeling, Mühlenstädt et al. (2012, p. 735) mention that the first case can lead to substantial modeling error whereas the second leads to a higher number of parameters, but still the true structure is contained in the model. In the optimization application, Ivanov and Kuhnt (2014, pp. 8-9) report that the first case of cutting an active interaction can 0 10000 20000 30000 40000 50000 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 number of evaluations TII Jansen Liu-Owen pick-freeze fixing method Figure 3.5: Simulation study on TII estimators, asymptotical exploration on the pure 3 rd order function. lead to a wrong optimum as the assumption of additivity is crucial in the parallelization whereas the keeping of an inactive interaction makes the procedure less efficient but does not effect the optimum. Finding at least all active interactions therefore is an important property of a good threshold.

In the following, a simple graphical method to support the threshold decision is introduced as well as a data-driven approach, which uses block-additive Kriging kernels.

Delta jump plot

As decision plots are a simple and clear way to help the user to find data dependent values, we suggest a plot, the delta jump plot, which helps to find suitable candidates The delta jump plot consists of two parts. In the bottom part, it shows the ordered values of all d 2 estimated TIIs. High jumps point to big differences between successive indices. In a situation with moderate perturbation, the difference between inactive and active indices should be high and thus be revealed by a jump. High jumps are additionally highlighted by shading of the relevant interval, the darker the shading the higher the jump, so that values of delta that lie inside dark intervals indicate good choices for the threshold. Simultaneously, the upper part of the plot shows the number of cliques of the corresponding FANOVA graph a threshold would lead to. A small number of cliques might be preferable in order to obtain a clearer structure. For the artificial case in Fig. 3.6, we see that a suitable δ lies between 0.0038 and 0.0196, for example 0.01, because the jump to the right of this interval is the highest and the number of cliques is very low. Another possible candidate would is 0.04.

Kriging kernel comparison

One method to compare different candidates for the threshold, possibly found by looking at the delta jump plot, is to compare the performance of Kriging models with corresponding block-additive structures.

As described in [START_REF] Mühlenstädt | Data-driven Kriging models based on FANOVA-decomposition[END_REF] and already mentioned in Section 3, Kriging kernels can be adapted to mimic the interaction structure given from the TII interaction screening. Instead of a product kernel (compare (2.3)), a block-additive kernel is applied. Basing on the cliques C 1 , . . . , C L of the FANOVA graph, it is constructed as the sum of kernels k C 1 , . . . , k C L , where each k C l is chosen to be a tensor-product kernel

in dimension |C l |, k(h) = L l=1 k C l (h C l ).
Kernels adapted to the true interaction structure of the underlying model are supposed to be more suited to represent the underlying model and thus lead to better predictions.

The idea for the threshold decision is to set up Kriging models according to some candidate thresholds and compare their leave-one-out prediction performance. See below for the precise approach.

• Choose a (low) number c of possible candidate thresholds and add 0 and 1 for the extreme cases of all input variables interacting and no interaction: δ = (0, δ 1 , . . . , δ c , 1).

• For each threshold δ i ∈ δ:

-Determine the according active interactions following rule (3.12).

-Set up the FANOVA graph and determine the cliques C 1 , . . . , C l .

-Construct the adapted Kriging model with block-additive kernel,

k(h) = L l=1 k C l (h C l ).
-Predict available original evaluations of the underlying model via leave-oneout cross-validations and compute the resulting RMSE value.

• Compare the leave-one-out RMSE values for each threshold candidate and choose the threshold with the smallest RMSE. It is also possible to plot scatterplots of predicted versus true values to get an impression of the precision of each prediction.

Simulation study

The performance of the threshold decision methods shall be presented in a small simulation study with four input variables (d = 4). Five different interaction structures (str1 to str5 ) are studied, here represented as sets of cliques of the corresponding graph structure. The goal is to find the threshold that rediscovers the corresponding graph.

• str1: {1}, {2}, {3}, {4} (pure additivity),

• str2: {1, 2}, {3}, {4} (two variables interacting),

• str3: {1, 2}, {3, 4} (two pairs of variables interacting),

• str4: {1, 2, 3}, {4} (three variables interacting, one additive),

• str5: {1, 2, 3, 4} (full graph).

For each structure, a number of 100 simulations are performed via the following steps, using mainly the R package fanovaGraph, described in the next chapter. A maximin distance design within the class of Latin hypercubes [START_REF] Morris | Exploratory designs for computational experiments[END_REF] is constructed for a number of 50 runs using the R package DiceDesign [START_REF] Franco | DiceDesign: Designs of Computer Experiments[END_REF]. As underlying model f , a Gaussian process with block-additive kernel corresponding to one of the interaction structures str1,. . . , str5 is simulated at the design points. See [START_REF] Fruth | The fanovagraph package: Visualization of interaction structures and construction of block-additive Kriging models[END_REF] for the simulation of block-additive processes. On the interaction structure str1 str2 str3 str4 str5 complete graph identified 0.99 0.68 0.50 0.85 0.60 all active interactions identified 1.00 0.98 0.81 0.95 0.60 Table 3.5: Proportions of correctly identified graphs in threshold simulation.

simulated data, an initial metamodel is fit. We use a Gaussian process model with a special ANOVA kernel [START_REF] Durrande | ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis[END_REF] since it is assumed to be more adapted to block-additive situations than the standard product kernel (2.3), which a priori assumes a full graph. For each interaction, the TII is estimated via the Liu and Owen TII estimator on a number of 20 000 evaluations for the whole set. The Kriging kernel comparison procedure described above is then executed. The candidate thresholds for the procedure are obtained from the three biggest jumps in the delta jump plot.

The results can be seen in Tab. 3.5. In the first row, the proportion of all correctly identified graphs is shown. The second row shows the proportion of those situations where all active interactions are identified correctly but not necessary all inactive ones.

As described before, this property is crucial for subsequent methods. A threshold that identifies all active interactions does at least not disturb subsequent analysis.

For str1 and str5, the results of the second row are obvious due to their interation structure. The simulation results show that the performance depends on the interaction structure. Situations with few interactions are better discovered than complex ones.

The interaction structure str3 with two similar interactions seems to be more difficult to predict than the unbalanced one str4. The results show that the Kriging kernel comparison leads to the right decision in the majority of the simulations and, except for the most complex case of a full graph, the active interactions are identified in more than 4 out of 5 simulations.

The simulations are repeated using the standard product kernel instead of the ANOVA kernel in the initial model in order to study the impact of the initial kernel on the interaction structure str1 str2 str3 str4 str5 complete graph identified 0.99 0.40 0.03 0.18 0.87 all active interactions identified 1.00 0.97 0.79 0.87 0.87 Table 3.6: Proportions of correctly identified graphs in threshold simulation using the standard product kernel as initial kernel.

performance. The results can be seen in Tab. 3.6. As assumed, the performance stays behind the performance of the ANOVA kernel, especially for the cases with two cliques, str3 and str4, for which only a very small proportion of the simulation is identified

correctly. An exception is the full graph str5, which was identified well. Despite the worse performance in the graph identification, the kernel lead to satisfactory results in the identification of all active interactions. This shows that with the standard kernel the procedure tends to identify too much interactions as active, explainable by the product structure of the kernel, which assumes a full graph. Thus, in this simulation study on simulated block-additive kernels, the ANOVA kernel performed generally better than the standard product kernel.

Implementation

As part of this work, the methods presented in this chapter are implemented in the R package fanovaGraph [START_REF] Fruth | The fanovagraph package: Visualization of interaction structures and construction of block-additive Kriging models[END_REF], which is online available on the repository CRAN (Comprehensive Archive Network). It includes the TII estimators of Section 3.2, the threshold decision methods of Section 3.6, and the FANOVA graph visualization mentioned in the motivational example together with Kriging model adaption. The package is explained in detail in the manual by [START_REF] Fruth | The fanovagraph package: Visualization of interaction structures and construction of block-additive Kriging models[END_REF].

Crossed DGSM

Similar to the total sensitivity index and the DGSM, the notion of the TII can as well be formulated as a derivative-based index, which provides a faster-to-compute index, which serves as an upper bound for the TII. In the context of statistical learning, [START_REF] Friedman | Predictive learning via rule ensembles[END_REF] introduced a quantity that generalizes the DGSM to interactions by using cross-partial derivatives and that we propose to denote by crossed DGSM,

ν i,j = ∂ 2 f (x) ∂x i ∂x j 2 dµ(x). (3.13) It is assumed that f fulfills ∂f (x) ∂x i ∂x j ∈ L 2 (µ).
From the definition, it is apparent that if ν i,j = 0, no term of the FANOVA decomposition (2.4) that contains x i and x j can be active at the same time, i.e. ν i,j = 0 ⇒ D i,j = 0. The crossed DGSM can thus be used for interaction screening.

Link to the TII

In [START_REF] Roustant | Crossed-derivative based sensitivity measures for interaction screening[END_REF] it is shown that the crossed DGSM provide upper bounds for the TII like the DGSM for the total sensitivity index.

Proposition 3.8. Let us consider n distributions µ 1 , . . . , µ n on the real line R, and

µ = µ 1 ⊗ • • • ⊗ µ n .
Assume that all µ i (i = 1, . . . , n) satisfy the Poincaré inequality (2.32). Let g : R n → R be a function in L 2 (µ), such that all first-order and crossed second-order partial derivatives are in L 2 (µ). Then for all pairs {i, j} (1 ≤ i, j ≤ n),

D i,j ≤ C(µ i )C(µ j )ν i,j . (3.14) Furthermore, C opt (µ i )C opt (µ j )

Estimation

As the crossed DGSM serve as upper boundaries for the TII, they can be used as a cheaper method for the screening of interactions, especially if the Hessian of the underlying model is available. As this is rarely the case, the necessary second-order derivatives usually have to be approximated, e.g. by finite differences. With x a n × d-data set and δ * a small real number, an estimator is given by

ν i,j = 1 n n k=1 ∆ ν (δ * ) 2 2 , with ∆ ν = f (x (k) i +δ * , x (k) j +δ * , x (k)
-{i,j} )-f (x

(k) i +δ * , x (k) j , x (k)
-{i,j} )-f (x

(k) i , x (k) j +δ * , x (k) -{i,j} )+f (x (k) ).
Analytical examples as well as practical applications can again be found in [START_REF] Roustant | Crossed-derivative based sensitivity measures for interaction screening[END_REF].

Conclusion

This chapter introduced a complete methodology for sensitivity analysis of the interaction structure of a black-box function as extension to a classical sensitivity analysis of first-order and total effects. The total interaction index (TII) was introduced, which, in the same way as the total sensitivity index screens out the most influential input variables, provides a screening of interactions. It reveals the way the input variables interact, which can be visualized in the FANOVA graph. The resulting information on the block-additive interaction structure can further be employed in adapted modeling or parallelized optimization.

Five possible estimators were presented together with their properties. The Liu and

Owen TII estimator was proven to have good properties to estimate a single TII, as it is unbiased, nonnegative, and asymptotically efficient. Its superiority was confirmed in simulations on analytical functions. Only for complex functions with a high number of inputs, the pick-freeze TII estimator performed better. There, the increase of costs when the dimension increases is linear instead of quadratic as for the Liu and Owen TII estimator. Furthermore, an approach to find a threshold for interaction screening was presented, a combination of a graphical decision and Kriging model comparison.

The corresponding R package fanovaGraph was presented. Finally, crossed DGSM were introduced as extensions of the DGSM to interaction analysis, and the inequality between total effects and DGSM was generalized in the context of interaction screening.

So far, all input variables were assumed to be continuous scalar variables. The next chapter presents a method to explore the sensitivity of additional functional input variables.

Chapter 4

Sensitivity Analysis for Functional Input

The complexity of the input settings of computer experiments is usually not limited to scalar numbers. Examples are industrial processes, where parameters can be varied over the process time or parameters in geological simulations, which have a spatial distribution. The inputs in these examples can be considered as functional input variables, dependend on time or space, respectively. They introduce an additional dimension to computer experiments and make new methods for the analysis necessary. For sensitivity analysis, it is especially interesting to explore the influence over the functional domain in order to find out about the inputs influence at different regions of the domain.

Usually, data with functional input is studied by functional linear regression (see Ramsay and Silverman (1997)), a framework that allows for approximation, modeling and prediction of data with functional input. Applied as a technique for sensitivity analysis, it has some drawbacks. In the functional linear regression framework, data is usually assumed as already provided so that no statistical design techniques are required. In addition, the interpretation of the influence of different regions of the functional domain is not easy (see e.g. James et al. ( 2009)) and it is restricted to linear behavior. Nonlinear modeling, however, is possible by functional Gaussian process frameworks as for instance presented by [START_REF] Morris | Gaussian surrogates for computer models with time-varying inputs and outputs[END_REF] and, specifically for scalar output, by [START_REF] Mühlenstädt | Computer experiments with functional inputs and scalar outputs by a norm-based approach[END_REF]. These methods also allow for functional space-filling design, but are not constructed for sensitivity analysis. [START_REF] Iooss | Global sensitivity analysis of computer models with functional inputs[END_REF], [START_REF] Lilburne | Sensitivity analysis of spatial models[END_REF]Taran-tola (2009), and[START_REF] Fort | Estimation of the Sobol indices in a linear functional multidimensional model[END_REF] eventually present several ideas for the sensitivity analysis of computer experiments with functional inputs. The methods result in one uncertainty index for each functional input as a whole and thus do not give insight into the sensitivity with respect to changes at specific intervals of the functional domain.

This chapter presents a general method for the performance of sensitivity analysis in computer experiments that can take functions as input variables, and return a scalar value as output. The aim is not only to discover the sensitivity of those functional variables as a whole but to identify relevant regions in the functional domain. A very economical sequential approach is presented, which reduces the functional space to a scalar one and makes use of methods from the so-called group factor screening. It results in a descriptive graphical representation of the functional sensitivities. The chapter starts with a motivational example. Then the methodology is presented, followed by properties and graphical presentation. Afterwards, the sequential design strategy for functional input is described. The chapter is finished with some advise on the practical implementation, including the R implementation as well as a detailed presentation of the procedure on an analytical example.

The results of the chapter are based on the contribution "Sequential designs for sensitivity analysis of functional inputs in computer experiments" by Fruth, Roustant, and Kuhnt (2014a).

Motivational example

To illustrate the problem, an application in sheet metal forming is introduced, which originally motivated this research as part of the PhD thesis. The complete analysis of the problem along with a detailed description follows later in Chapter 6.2. Here, it suffices to say that the aim is to analyze the sensitivity of a scalar output, springback, for two inputs that can be changed during the forming process: friction and blankholder force. The corresponding computer experiment is very time-consuming, which is even enhanced through the functional variation of the two parameters.

A few questions arise at this point. Does the effort of varying blankholder force and friction during the process indeed change anything in the outcome? For example, might it be that only the area under each input curve, and thus the mean of the settings over the time, has an influence? In this case, the functional approach would be useless as it would not give more information than the scalar analysis, where each input variable is kept constant during the run.

To check this issue, a small preview study on a few simple runs is executed (Fig. 4.1).

Each row shows the setting of the parameter friction of one run of the computer experiment. No other parameter is changed over the runs, so that the functional behavior of the friction can be analyzed exclusively. In the first two runs, the friction is constant over the 15 seconds of the punch travel, leading to a springback of 7.1 for a low and 2.2 for a high value. In the following 4 runs, varying curves are chosen such that the area below each curve is the same in each run. It can be seen that varying the friction can dramatically reduce springback. In addition, although the four last runs have equal mean friction over time, very clear differences in the springback results are obtained. This indicates the potential that lies in the functional approach and motivates an exploration of the function influence of the parameters.

Notation and situation

For the handling of functional inputs, the notation is slightly modified. An underlying model f is considered with d scal scalar input variables functional inputs. There are no further conditions on the functions, that is they can be designed freely with no assumptions on the shape or the smoothness of the functions.

x i ∈ [-1, 1], i = 1, . . . ,
The input parameters are connected to Y by f via f :

[-1, 1] d scal × F d fun [0,1] → R, Y = f (x 1 , . . . , x d scal , g 1 , . . . , g d fun ) ,
where F [0,1] denotes the space of all functions on [0, 1].

Sequential functional sensitivity analysis

Modeling of input functions

With functional inputs, the input dimension becomes infinite dimensional while at the same time the computer experiment becomes more time-consuming. The functions shall therefore be constructed within a flexible framework that reduces the input space considerably, and that eventually allows for the identification of important regions of the functional domain. The idea is to use piecewise constant functions as shown in Fig. 4.2, that is to restrict the input functions to the space of piecewise constant functions. Say we have a decomposition of the domain of each functional input g j into p j ∈ N + subintervals at splitting points a j = (a 0 j , . . . , a

t a 0 = 0 a 1 a 2 a 3 a p-2 a p-1 a p = 1 (g i ) t -1 1 Z (1) i Z (2) i Z (3) i Z (p-1) i Z (p) i . . . . . .
p j j ) with 0 = a 0 j < a 1 j < • • • < a p j -1 j < a p j j = 1, D j = [a 0 j , a 1 j [ ∪ [a 1 j , a 2 j [ ∪ . . . ∪ [a p j -1 j , a p j j ].
We restrict each g j to belong to V a j , the space of piecewise constant functions over the subintervals defined by a j , V a j = Z

(1)

j 1 [0,a 1 j [ (t) + • • • + Z (p j ) j 1 [a p j -1 j ,1] (t), with Z (k) j ∈ [-1, 1] , 1 ≤ k ≤ p j . (4.1)
The specific procedure of how to choose the split points is provided later in Section 4.3.

With this representation, every functional input g j ∈ V a j is identified by the random variables Z

(1) j , . . . , Z (p j ) j

. Hence, the input dimension is transformed from functional to scalar space, Y =f (x 1 , . . . , x d scal , g 1 , . . . , g d fun )

= fa 1 ,...,a .

The piecewise constant representation can be viewed in the context of functional representation, see e.g. [START_REF] Ramsay | Functional data analysis[END_REF], Chapter 3. B-Splines, piecewise polynomial functions with certain knots as breakpoints, are a popular basis to model functions [START_REF] De Boor | A practical guide to splines[END_REF]. They cover various types of functions and are always bounded. Our functional framework (4.1) can be regarded as a B-spline of order 1, a linear combination of piecewise constant functions. In addition, the transformation is connected to wavelet theory, more precisely to Haar Wavelets (see e.g. [START_REF] Walker | A primer on wavelets and their scientific applications[END_REF]),

where again the basis functions are piecewise constant functions but scaled and shifted to the desired functional form. Wavelets and our framework resemble each other also in the sequential splitting procedure described later in Section 4.3, where sequences of embedded spaces V a j ⊆ V a j ⊆ V a j ⊆ . . . are considered. However, the objectives of this sequential splitting are different, as in our approach, the interest lies in the time localization, the identification of influential time points, not in the time-frequency localization.

Remark: Scalar inputs The design of scalar inputs does not have to be mentioned separately, as the space of piecewise constant functions also includes scalar inputs. A scalar input can be considered as a functional input over only one interval [0, 1] which is never split. Thus, from now on, d = d scal + d fun is simply used to denote the number of inputs.

Sensitivity Indices

The defined framework allows to perform sensitivity analysis on the values

Z (k) j in Y = fa 1 ,...,a p d Z (1) 1 , . . . , Z (p 1 ) 1 , . . . , Z (1) d , . . . , Z (p d ) d
, leading to a clear and easily understandable sensitivity analysis of intervals. The specific sensitivity method can be chosen according to the evaluation costs and complexity of the problem from the methods described in Chapter 2. However, since the functional computer experiment is assumed to be expensive, regression analysis as a very economical method for sensitivity analysis is recommended and explored further in this work.

The following linear model is considered, Additionally, higher order effects of interest can be estimated. One might be especially interested in the second-order interactions between intervals of the same functional input, which leads to the model

Y = α + d j=1 p j k=1 β (k) j Z (k) j + ε, ( 4
Y = α + d j=1 p j k=1 β (k) j Z (k) j + d j=1 1≤k<k ≤p j β (k,k ) j Z (k) j Z (k ) j + ε. (4.3)

Graphical representation and normalization

The functional indices can be represented by plotting them over the functional domain of each functional input, which visualizes the functional behavior of each input in a compact way. An artificial example of such a visualization is depicted in Fig. 4.3 (top left). It shows the regression coefficients of a single functional input, the bar colors highlight direction and size of the index. A positive bar size indicates that a raise of the input at this interval causes an increase in the mean output, a negative size causes a decrease. In the specific plot, a small negative influence is visible for the first half of the functional domain and a larger positive influence in the second half. In a study with more than one functional input, one such plot is drawn for each input. Scalar inputs can be depicted as single bars.

The drawback of using β directly as functional sensitivity index is that the index values are not independent of the current splitting a j = (a 0 j , . . . , a p j j ). If we change the splitting so that the last two intervals are gathered, this larger interval now has a much stronger influence than the two single intervals before, since a larger part of the functional input is varied, which leads to a larger change in the output. is defined by

H k j = β (k) j a k j -a k-1 j ,
and the normalized interaction regression index of

Z (k) j and Z (k ) j by H k,k j = β (k,k ) j (a k j -a k-1 j )(a k j -a k -1 j ) for j ∈ {1, . . . , d}, 1 ≤ k < k ≤ p j .
In picture on the top. In addition, the bar of the first interval is now smaller, as its size is larger and thus was not comparable to the other intervals in the upper plots.

As a regression coefficient β can be interpreted as the change in the output when the functional input is increased over the corresponding interval, H can be interpreted as the change in the output when increasing the functional input over one unit of the functional domain, e.g. one time unit for temporal inputs. Yet, care has to be taken in the interpretation. In an extreme case, the joint interval of two intervals with influences of different signs could show zero influence when the intervals cancel each other out.

Any interpretation of the indices can only relate to the mean influence of the interval, as can be seen in the following test cases.

Error-free test cases

In order to show that the normalized indices fulfill desirable properties, their performance on some error-free test cases for the underlying model f shall be examined. Each test case depends on only one functional input g and follows a purely linear behavior.

Proposition 4.1. Let f be the integral over the functional input g weighted by an

integrable function w : [0, 1] → R f (g) = α + 1 0 w(t)g(t) dt,
with α ∈ R. Then, for any splitting a = (a 0 , . . . , a p ) of the functional domain of g, the normalized regression indices return the mean values of the weight function over the intervals, and each interaction index is zero,

H k = a k a k-1 w(t) dt a k -a k-1 and H k,k = 0, k, k = 1, . . . , p.
Asymptotically, if the interval size goes to zero, the normalized regression indices return the exact value of the weight function,

lim ∆→0 t+∆ t w(u) du ∆ = w(t).
Proof. As g is bounded and w is integrable, wg is integrable. When we insert the piecewise constant function representation

g(t) = p k=1 Z (k) 1 [a k-1 ,a k [ (t)
into f (g), we obtain

f (g) =α + 1 0 w(t) p k=1 Z (k) 1 [a k-1 ,a k [ (t) dt =α + Z (1) a 1 a 0 w(t) dt + • • • + Z (p) a p a p-1 w(t) dt =(1, Z (1) , . . . , Z (p) )η, with η = 1, a 1 a 0 w(t) dt, . . . , a p a p-1
w(t) dt . Now in the experiments, recalling that the design matrix Z on the values of Z (1) , . . . , Z (p) is such that (Z Z) is invertible, and denoting by y the corresponding output, we get y = Zη. The Least Squares coefficients are then given by

β = (Z Z) -1 Z y = η,
which leads to

H k = a k a k-1 w(t) dt a k -a k-1 .
As there are no interaction terms in Z, we get H k,k = 0, k, k = 1, . . . , p..

In the same way it can be shown that interactions are recovered.

Proposition 4.2. Let f be the product of two integrals of g over two different intervals

[a i * -1 , a i * ] and [a j * -1 , a j * ] with c, α ∈ R f (g) = α + c a i * a i * -1 g(t) dt × a j * a j * -1 g(t) dt.
Then for any splitting a of the functional domain of g with a i * -1 , a i * , a j * -1 , a j * ∈ a, the estimates of the indices and the interaction indices recapture the interactions as they hold

H k = 0, k = 1, . . . , p and H k,k =      c, k = i * , k = j * , 0, otherwise.
Proof. Similarly to the proof of Prop. 4.1 we get

f (g) =α + c a i * a i * -1 g(t) dt × a j * a j * -1 g(t) dt =α + c a i * a i * -1 p k=1 Z (k) 1 [a k-1 ,a k [ (t) dt × a j * a j * -1 p k=1 Z (k) 1 [a k-1 ,a k [ (t) dt =α + c × Z (i * ) (a i * -a i * -1 ) × Z (j * ) (a j * -a j * -1 ).
Then the Least Squares estimation returns β k = 0, k = 1, . . . , p and

β k,k =      c(a i * -a i * -1 )(a j * -a j * -1 ), k = i * , k = j * ,

0, otherwise

leading to

H k = 0, k = 1, . . . , p and H k,k =      c, k = i * , k = j * , 0, otherwise.
Furthermore it can be shown that the normalized regression index is robust against nonlinear transformations of the input g in the computer model.

Proposition 4.3. If g is transformed by a not necessarily linear, but strictly mono-

tonically increasing function ζ : [0, 1] → R f (g) = α + 1 0 w(t)ζ(g(t)) dt
then for a two-level design the values H k are linear in the mean values of the weight function, and the weighting function still determines the importance of the intervals compared to each other. Thus, there exists a fix λ > 0 such that

H k = λ a k a k-1 w(t) dt a k -a k-1 , k = 1, . . . , p.
Proof.

f (g) =α + 1 0 w(t)ζ(g(t)) dt, =α + ζ(Z (1) ) a 1 a 0 w(t) dt + • • • + ζ(Z (p) ) a p a p-1 w(t) dt.
Now for a design on two levels {-1,1}, we can define λ and κ ∈ R, such that

ζ(z) = λz + κ for the two values z = -1, z = 1.
Since ζ is increasing, we have λ > 0. It follows for all

Z (k) ∈ {-1, 1}, k = 1, . . . , p that f (g) =α + (λZ (1) + κ) a 1 a 0 w(t) dt + • • • + (λZ (p) + κ) a p a p-1 w(t) dt =   α + κ p k=1 a k a k-1 w(t) dt   + p k=1   λ a k a k-1 w(t) dt   Z (k) .
Then in the computations of the coefficients, we obtain

β k =λ a k a k-1 w(t) dt, k = 1, . . . , p ⇒ H k =λ a k a k-1 w(t) dt a k -a k-1 , k = 1, . . . , p.
We have now described the functional shape for the design of the input functions, the way to perform sensitivity analysis, and its visualization. In the following, the missing part on how to choose the splitting into the intervals and appropriate design techniques are presented.

Splitting and design

If a desired splitting is available, the method can be performed right away, using a suitable design matrix Z, e.g. full or fractional factorial design. If not, one aims at discovering the functional domain while investing as few as possible model evaluations. Therefore an economical splitting approach is proposed which is performed in consecutive steps. In each step preceding information is used for the splitting.

For the sake of readability, only one functional input g : D → [-1, 1] as sole input is considered in this section, i.e. d = 1. The approach is easily extended to more inputs by considering them as additional groups. For a specific iteration step r, we slightly change the notation of the splitting points to a r = (a 0,r , . . . , a p r ,r ), 0 = a 0,r < a 1,r < • • • < a p r -1,r < a p r ,r = 1 the space of piecewise constant functions to

V a r = Z (1,r) 1 [0,a 1,r [ (t) + • • • + Z (p r ,r) 1 [a p r -1,r ,1] (t), Z (k,r) ∈ [0, 1] ,
and the normalized indices to

H k r , k = 1, . . . , p k and H k,k r j , 1 ≤ k ≤ k ≤ p k .

Sequential splitting approach

The approach is based on a family of methods called group factor screening (see [START_REF] Watson | A study of the group screening method[END_REF] or for an overview Morris ( 2006)), a very economical method to screen influential input variables in experiments with a high number of input variables. The basic idea is to group variables, explore the influence of the groups as a whole, and sequentially divide only those groups that are influential. In the literature on group factor screening, there are various ways on how to design the groups at the different steps, which differ e.g. in their assumptions, orthogonality, the treatment of interactions, or the way of reusing evaluations.

The idea can be transferred to functional sensitivity analysis by interpreting the (infinite) points in the functional domain as individual input variables and the splitting intervals as groups of them. The approach is then to start with a very low number p 1 of intervals a 1 = (a 0,1 , . . . , a p 1 ,1 ), perform sensitivity analysis by choosing a suitable design for the corresponding variables Z (1,1) , Z (2,1) , . . . , Z (p 1 ,1) , perform the experiments, and compute the corresponding indices H 1 1 , . . . , H p 1 1 . Then, for the second step, only those intervals that show a considerable influence on the output are examined further. So the vector of splitting points of the second step a 2 includes all points of a 1 plus the point a k-1,1 +a k,1 2 for each interval a k-1,1 , a k,1 to be examined. This procedure of estimation and splitting is repeated until the functional domain is sufficiently explored or the maximum budget is reached. The split decisions should be taken in close cooperation with experts. Generally, an interval should be seen as important if its index is bigger than an assumed approximation error.

Remark: Interpretation care It has to be kept in mind that the indices are only mean values over a given interval, as it could be seen in Prop. 4.1. Intervals, even if the sensitivity index is small, should be explored further if a change in sign is suspected.

Design by sequential bifurcation

The application of a specific group factor screening method, the sequential bifurcation [START_REF] Bettonvil | Factor screening by sequential bifurcation[END_REF], shall be presented in detail here. In this method, the evaluations of the different steps are effectively reused in the subsequent steps, resulting in a very economical procedure.

Adapted to functional sensitivity analysis, the procedure is the following. The variables are designed on two extreme levels, encoded with -1 and 1. It starts with only one interval, corresponding to a constant curve, which is set to -1 and, in a second run, to +1. Using the results of these two runs, y 1 and y 2 , the normalized regression index of the whole function can be estimated by H 1 1 = 0.5(y 2 -y 1 ) 1 . The domain is then split in the middle, resulting in two intervals in the second step. Only one additional run is required to estimate the coefficients of both intervals in which the first interval is set to +1 and the second interval to -1, resulting in y 3 . Values for (+1, +1), (-1, -1) are already known from the previous step, so that the indices of both intervals can be estimated by H 1 2 = 0.5(y 3 -y 1 ) 0.5 and H 2 2 = 0.5(y 2 -y 3 ) 0.5

t 0 1 Z (1,1) = -1 run 1 → y 1 Z (1,1) = +1 run 2 → y 2 Z (1,2) = +1 Z (2,2) = -1 run 3 → y 3 Z (1,3) = +1 Z (2,3) = -1 Z (3,3) = -1 Z (4,3) = -1 run 4 → y 4 Z (1,3) = +1 Z (2,3) = +1 Z (3,3) = +1 Z (4,3) = -1 run 5 → y 5 . . . . . . . . .
. Generally, for each split, only one additional run is required to estimate the influence of both new intervals, the one with +1 up to the cut and -1 from there. The approach is depicted in Fig. 4.4. It shows three sequential steps requiring a total of 5 runs.

When the presence of interactions has to be considered, additional mirror runs are suggested by [START_REF] Bettonvil | Factor screening by sequential bifurcation[END_REF], i.e. adding a new run for each run in the design that contains the same settings but with opposite signs. By this, unbiased estimates of the first-order effect coefficients are obtained.

Different designs may be necessary if orthogonality and/or the possibility to estimate interactions are required. Then factorial or fractional factorial designs are good design options. A new such design is then constructed in each step r on all variables of interest Z (1,r) , . . . , Z (p r ,r) , all noninteresting intervals are set to a constant value, e.g. 0, in the design and thus are not regarded in the sensitivity analysis any more. Here again, runs from former steps can be reused. It is easy to show that for full factorial designs in a step where all intervals have been split half of the required runs can be reused. For fractional factorial designs, the reuse possibilities depend strongly on the confounding.

Implementation

The presented methodology is implemented in the R package seqSAFI [START_REF] Fruth | seqSAFI: sequential sensitivity analysis for functional inputs[END_REF]. The package allows for all described steps of sequential design, modeling and plotting of the functional sensitivities. See

Analytical example

To show how the method can practically be applied, a small analytical test example is analyzed in this section. For the unknown and expensive to evaluate underlying model, we choose the following function which depends on two functional inputs, g 1 and g 2 ,

f (g 1 , g 2 , x) = 1 0 (3 -9t)1 [0, 1 3 ] (t) + Z(t) g 1 (t) dt - 1 0 1 30 1 [ 3 10 ,1] (t)g 2 (t) + 3 3 dt + 8 10 sin(x), (4.4)
with Z(.) a simulation of a Gaussian process with zero mean, a small variance of 0.05, and Matérn 5/2 covariance function with parameter θ = 0.1. The process Z(.) as well as the third power and the sine function are disturbances meant to achieve a more realistic example. The aim is to discover the linear effects of each input, that is

• the decreasing positive influence of g 1 over 0, 1 3 , • the constant negative influence of g 2 over 3 10 , 1 , • the positive influence of x.

The computations are performed using the R package seqSAFI, presented before. At first, a starting design has to be created. As we want to use as little evaluations as possible and we say that we do not assume active interactions, we use a sequential bifurcation design without mirror runs. Not having any previous knowledge about the functional influence, we start with two equidistant intervals per functional input, i.e with the decompositions

a 1 = a 2 = 0, 1 2 , 1
for g 1 and g 2 respectively. The scalar input x is treated as a functional input on the sole interval a 3 = (0, 1) in all steps. Thus, a starting design is set up that contains six runs, the initial bifurcation runs (all +1 and all -1) plus one run for each splitting which makes two for the splitting into three variables and two for the splitting of each of the two functional domains. The design can be seen on the top in The intervals with visible influence of the second step are split again to gain a more detailed look in the functional influence, leading to two new runs for g 1 and three for g 2 . They are shown in Fig. 4.8 together with the resulting indices. We stop the procedure at this point, making Fig. 4.8 the final result. Though they are disturbed and do not show the exact influences, the main linear behavior of the three parameters is revealed in the 15 runs. We can see the decreasing influence of g 1 at the beginning of the domain, though not exactly at the first third, the constant influence of g 2 in the last part, and the positive influence of the scalar input x.

input g 1 g 2 x interval 0, 1 2 1 2 , 1 0, 1 2 1 2 , 1 ([0, 1]) run 1 +1 +1 +1 +1 +1 run 2 -1 -1 -1 -1 -1 run 3 +1 +1 -1 -1 -1 run 4 +1 +1 +1 +1 -1 run 5 +1 -1 -1 -1 -1 run 6 +1 +1 +1 -1 -1 g 1 normalized regression index

Chapter discussion

In this chapter a methodology was presented to perform sensitivity analysis for timeconsuming computer codes that take, additional to scalar inputs, functions as input variables. The focus was not only on discovering the sensitivity of the function as a whole but also on exploring the behaviour over the functional domain. The method comprises sequential group factor screening of piecewise constant functions and specially normalized regression indices.

Due to the infinite functional dimension and the small budget of evaluations, a very economical method has been developed that has several limitations. Independence as well as only linear behaviour over the intervals is assumed. Another limitation concerns the interpretation of the indices, as theoretically influences of opposite sign could get canceled out. The method, however, was shown to be robust against some nonlinear transformations and proves to work well in the application in Chapter 6.

input g 1 g 2 x interval 0, 1 6 1 6 , 1 3 1 3 , 1 2 1 2 , 1 0, 1 4 1 4 , 1 2 1 2 , 3 4 , 1 ([0, 1]) run 7 +1 -1 -1 -1 -1 -1 -1 -1 -1 run 8 +1 +1 -1 -1 -1 -1 -1 -1 -1 run 9 +1 +1 +1 +1 +1 -1 -1 -1 -1 run 10 +1 +1 +1 +1 +1 +1 +1 -1 -1 g 1 normalized regression index
run 11 +1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 run 12 +1 +1 +1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 run 13 +1 +1 +1 +1 +1 +1 +1 +1 -1 -1 -1 -1 -1 -1 run 14 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 -1 -1 -1 -1 run 15 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 +1 -1 -1 g 1 normalized regression index
The basic idea of this chapter, to explore single intervals and reduce their size more and more, can be transferred to another application of sensitivity analysis, the analysis of the support of input variables, which is presented in the next chapter.

Chapter 5 The first and the last variable have clearly lost influence in comparison to the second one. This can be explained by the interaction between X 1 and X 3 , which is much stronger at the borders than in the rest of the function. The index of X 2 stays unchanged as it does not interact with X 3 . The change in the output when changing X 2 stays the same.

Support Analysis

As very cheap-to-evaluate metamodels are often used in practice (Section 2.2), it is possible to compute more than the common one or two simple numbers, the Sobol and the total sensitivity index, for each input variable. In this chapter, an index that extends the traditional global sensitivity analysis is defined and examined. It explores the influence of different parts of the support of the input variables by applying the ideas of functional sensitivity analysis from Chapter 4 to a regular scalar situation. In a way, this method is a turn back to local sensitivity analysis (see Section 2.1), where, out of necessity, only the local influence of parameters is studied. The main difference is that in local sensitivity analysis the influence of a variable is studied in a local space of the other inputs whereas the following method studies the local influence of a variable, but in a global space of the other variables. In addition, the local influences are studied over the whole input space, so no space is ignored. This idea, more specifically the idea to use sensitivity analysis to determine regions in the input space for which the model variation is maximum, is mentioned in Saltelli et al. (2000, p. 6), but to our knowledge not pursued further.

The following method presents a way to show the behavior of the sensitivity for different parts of the input support, so that the effects of the choice of an input distribution can be visible. In addition, it gives further insight into the function by presenting a function of sensitivities as extension to Sobol and total sensitivity indices. Following the definition of the support index functions, we show that they are connected to classical indices when the variable in question is restricted to an interval whose size goes to zero. Then, the method is applied to the Ishigami function, followed by a study of their expected value over the support. The support index functions are then defined as follows.

Support index functions

Definition 5.1. Support index functions

The first-order support index D i (t) of an input variable X i at a point t, t ∈ [0, 1], is defined as the square of the expected value of the first derivative of f ,

D i (t) = E ∂f ∂x i (t, X -i ) 2 .
The total support index D T i (t) of an input variable X i at a point t, t ∈ [0, 1], is defined as the expected value of the squared derivative of f ,

D T i (t) = E ∂f ∂x i (t, X -i ) 2 .
The support variance D X i (t) corresponding to a support index of an input variable X i at a point t, t ∈ [0, 1], is defined as the overall variance for

X i = t, D X i (t) = Var(f (t, X -i )).
The functions can be evaluated and plotted for a sufficiently large number of discrete points t = 0, . . . , 1 over the support. If the gradient of f is known the functions can be calculated directly. If not, the gradient can be approximated by finite differences.

For instance the first-order support index of X i at point t using Monte Carlo samples of X -i , x

-i , . . . , x

-i , and a small value of δ * can be estimated by

D i (t) = 1 n n k=1 f (t + δ * , x (k) -i ) -f (t -δ * , x (k) -i ) 2δ * 2 .
An interpretation of the support functions can be obtained, when we look at the firstorder and total influence of the specific point in the support of the variable. Indeed, the functions are connected to the classical index estimation when we restrict X i to vary only over a small interval around t and let the size of the interval go to zero.

Proposition 5.1. Let X 1 , . . . , X d be independent random variables with continuous density functions with support [0, 1]. Let f be a function defined on ∆ = [0, 1] d compact and assume that f is of class C 2 , ensuring that f (X), ∂f ∂x• (X), ∂ 2 f ∂x•∂x• (X) are in L 2 . Now let one of the random variables, say X i , be restricted to

X h i ∼ X i |X i ∈ [t -h /2, t + h /2] for t ∈]0, 1[ and h > 0.
Then we have for the overall variance, the Sobol index of X h i , and the total sensitivity index of X h i that

1. Var(f (X h i , X -i )) h→0 --→ Var(f (t, X -i )) = D X i (t), 2. Var E f (X h i , X -i )|X h i /h 2 h→0 --→ 1 12 E ∂f ∂x i (t, X -i ) 2 = D i (t)/12, 3. E Var f (X h i , X -i )|X -i /h 2 h→0 --→ 1 12 E ∂f ∂x i (t, X -i ) 2 = D T i (t)/12.
Proof. Let us first show that

V ar(X h i -t) h 2 h→0 --→ 1 12 .
Let g be the density function of

X i -t, g(x) = f X i (x + t).
By assumption, g is continuous in 0, g(0) > 0. Consider the density of the truncated variable

X h i -t, g h (x) = g(x) h 2 -h 2 g(t) dt 1 [-h 2 , h 2 ] (x). Let h = inf x∈[-h 2 , h 2 ] g(x) and u h = sup x∈[-h 2 , h 2 ] g(x)
. With continuity of f in 0, we have

u h h→0
--→ g(0) and h h→0 --→ g(0). By standard integral computations we get

h h ≤ h 2 -h 2 g(x) dx ≤ hu h and h 3 12 h ≤ h 2 -h 2 g(x)x 2 dx ≤ h 3 12 u h so that for E X h i -t 2 = h 2 -h 2 g(x)x 2 dx h 2 -h 2 g(x) dx it holds that h u h h 3 /12 h ≤ E X h i -t 2 ≤ u h h h 3 /12 h .
By dividing by h 2 this leads to

E X h i -t 2 h 2 h→0 --→ 1 12 .
It remains to show that

E(X h i -t) h h→0 --→ 0. As h 2 -h 2 g(x)x dx = h 2 0 (g(x) -g(-x)) x dx it holds that |E(X h i -t)| ≤ 1 h h h 2 0 |g(x) -g(-x)| x dx ≤ 1 h h h 2 0 (u h -h ) x dx = u h -h h h h 2 0 x dx = u h -h h h h 2 8 . Thus |E(X h i -t)| h ≤ u h -h 8 h h→0 --→ 0.
In the following, denote B i =

X h i -t h . Remark that B i is bounded by |B i | ≤ 1 2 and that Var(B i ) h→0 --→ 1 12 .
Proof of 1. Write the mean value theorem between a real number

x i ∈ [t-h/2, t+h/2] ⊆ [0, 1] and t f (x) = f (x i , x -i ) = f (t, x -i ) + (x i -t) ∂f ∂x i (c, x -i ) for a c between x i and t, thus c ∈ [t -h/2, t + h/2].
Replacing by random variables, we get

f (X h i , X -i ) = f (t, X -i ) + (X h i -t) ∂f ∂x i (C, X -i ),
where

C is a random variable such that C ∈ [t -h/2, t + h/2]. Denoting Q = B i ∂f ∂x i (C, X -i ), we then have f (X h i , X -i ) = f (t, X -i ) + hQ.
Now remark that Q is bounded since both B i and ∂f ∂x i are bounded (by continuity of

∂f ∂x i on the compact set [0, 1] d ). This implies that E(f (X h i , X -i )) h→0 --→ E(f (t, X -i ))
by the Lebesgue dominated convergence theorem as well as E f (X h i ,

X -i ) 2 h→0 --→ E (f (t, X -i ) 2
). The result follows.

Proof of 2. The proof is similar to the one of 1. Write the Taylor-Lagrange expansion of f between a real number

x i ∈ [t -h/2, t + h/2] ⊆ [0, 1] and t f (x) = f (x i , x -i ) = f (t, x -i ) + (x i -t) ∂f ∂x i (t, x -i ) + 1 2 (x i -t) 2 ∂ 2 f ∂x 2 i (c, x -i )
for a c between x i and t, thus c ∈ [t -h/2, t + h/2]. Replacing by random variables, we have

f (X h i , X -i ) = f (t, X -i ) + hB i ∂f ∂x i (t, X -i ) + h 2 R, with R = 1 2 B 2 i ∂ 2 f ∂x 2 i (C, X -i ), where C is a random variable, C ∈ [t -h/2, t + h/2].
Then, using the independence between X h i and X -i , we have Proof of 3. With the same notations as in the proof of 2., we have

E f (X h i , X -i )|X h i = β 0 + hβ 1 B i + h 2 R, with β 0 = E (f (t, X -i )), β 1 = E ∂f ∂x i (t, X -i ) , and R = E R|X h i . Finally Var E f (X h i , X -i )|X h i = h 2 β 2 1 Var(B i ) + 2h 3 β 1 Cov(B i , R) + h 4 Var(R).
E f (X h i , X -i )|X -i = f (t, X -i ) + hE(B i ) ∂f ∂x i (t, X -i ) + h 2 E R|X -i and thus f (X h i , X -i ) -E f (X h i , X -i )|X -i = h(B i -E(B i )) ∂f ∂x i (t, X -i ) + h 2 S,
with S = R -E R|X -i . Since S is centered, we have by independence of B i and

∂f ∂x i (t, X -i ) E f (X h i , X -i ) -E f (X h i , X -i )|X -i 2 = β 2 Var(B i )h 2 + 2h 3 Cov (B i -E(B i )) ∂f ∂x i (t, X -i ), S + h 4 Var(S), with β 2 = E ∂f ∂x i (t, X -i ) 2
. The result follows in the same way as in the proof of 2. by using that ∂f ∂x i (t, X -i ) and S are bounded random variables.

Analytical example

Applying the described support analysis to the Ishigami function, sin(X 1 )+7 sin 2 (X 2 )+ 0.1X 4 3 sin(X 1 ), reveals the effect observed in the introduction. For the estimation, the finite element estimators D i (t) and D T i (t) were computed at 40 regularly spaced points in ]-π, π[ for Monte Carlo samples of size n = 5 000. See Fig. 5.1 for the result. The influence of the input variable X 3 is increasing strongly at the borders, which explains the observed strong change in the total effect, when the support size is changed. Beside this the sinusoidal behavior of X 1 and X 2 can be observed as well as the facts already known from the standard analysis that X 1 and X 3 interact and that the main influence of X 2 is the strongest.

Expected values

Taking the expected value over the support index functions reveals interesting connections to other indices, which give a little more insight into the interpretation of the functions.

For the total support index D T i (t), the expected value is equal to the DGSM ν

i (Sec- tion 2.4.2), E D T i (X i ) = E E ∂f ∂x i (X i , X -i ) 2 X i = ν i .
It follows that the expected value of the total support index serves as an upper bound for the total sensitivity index,

C(µ i )E D T i (X i ) ≥ D T i ,
with C(µ i ) the Poincaré constant of µ i as before.

∂f ∂x i 0 = E ∂f ∂x i (X) .
Now compute the first-order Sobol index of ∂f ∂x i (.) corresponding to the input variable x i by using the pick-freeze formula,

D i ∂f ∂x i (X) =E ∂f ∂x i (X i , X -i ) ∂f ∂x i (X i , Z -i ) - ∂f ∂x i 2 0 =E E ∂f ∂x i (X i , X -i ) X i 2 - ∂f ∂x i 2 0 . Thus we obtain E [D i (X i )] = D i ∂f ∂x i (X) + ∂f ∂x i 2 0 ,
the expected value over the first-order support index function can be obtained as firstorder index of the derivative function plus its global mean.

Chapter discussion

Two new indices, the first-order support index and the total support index, have been presented that extend the Sobol index and the total sensitivity index, respectively.

They return a function of sensitivities that give insight into the local behavior of input variables, and can be especially helpful in the specification of the input distribution.

It was shown that the support index functions are connected to the limits of their corresponding scalar indices when the support approaches zero. A link to DGSM and variance-based indices was made.

Chapter 6

6. Application to a Sheet Metal Forming

Process

The industrial process of deep drawing is a fundamental procedure for forming sheet metal into desired shapes. It is for example extensively applied in the automotive industry in the production of car bodies. In the deep drawing process, a flat sheet metal is pressed with a punch into a die while so-called blankholders keep the metal fixed at the metal borders. Three time points of the forming process are illustrated in Fig. 6.2 on the left. Examples of the target workpiece, inspired by the shape of B-pillar, can be seen in Fig. 6.2 on the right. One of the main problems to investigate in the deep drawing analysis is the shape accuracy of the formed part. Tearing or wrinkling can occur during the forming process. Another serious problem is springback, the reforming of the flange after the removal of the punch, which can lead to strong difficulties in the assembly. Such inaccuracy is exemplarily depicted in Fig. 6.2 on the right, springback can be observed in the right piece. The task of analyzing and optimizing the shape accuracy is growing in importance and complexity with the increasing requirements on the formed parts, especially in the automotive industry. In addition, advanced high strength steels are used more and more due to their advantages of light weight and higher strength, which are even more prone to springback.

The forming process can be adjusted by several parameters, which concern for instance the specific properties of the formed material, like thickness or strength, parameters that can be adjusted during the process like the blankholder force or the speed of the punch, or the geometry of the die. It is common practice to adjust and study the process on Finite-Element models, for instance by the popular Software tool AutoForm [START_REF] Autoform | AutoForm users manual[END_REF]. In this work, the more specialized Software LS-DYNA (Livermore Software Technology Corporation, 2005) is used for the deep drawing simulations. All other computations in the study are performed by the software R (R Core Team, 2014).

The aim of the study is to show how the forming process can be analyzed by the methods presented in this work, how they help to gain a deeper insight into the process and improve subsequent analysis techniques.

Thickness reduction

The first application aims at analyzing the influence of several input parameters on the thickness reduction of the formed part. A high thickness reduction at a region of the sheet metal means a strong local thinning, which can lead to structural deformations or actual tearing. Is it thus of interest to know, which parameters influence the thinning and in which way as well as modeling and optimization. The part under study is the demonstrator part already shown in Fig. 6.2 on the right. Eight input variables are varied in this study, listed in Tab. 6.1 together with their ranges. The variables flow stress, initial sheet thickness, hardening exponent and sheet layout concern the material, blankholder force and friction can be varied during the process. In practice, friction can be changed during the process by adding or removing lubricant with high pressure air and oil removing agents. To roughly integrate this change of friction over the process time, friction is modeled as three independent variables corresponding to three stages of the process time over which friction is kept constant. The more elaborate approach to the analysis of functional inputs of Chapter 4 is applied later in Section 6.2.

The output value is the maximal thickness reduction, that is the thickness reduction of the point of the formed part with the strongest thinning. It is given as the ratio between thickness reduction and initial thickness which leads to a scalar value between 0 and 1.

input variable feasible region

X 1 flow stress 100-200 MPa X 2 initial sheet thickness 0.5-1.7 mm X 3 blankholder force 50-200 kN X 4 friction; 1 st third of process time 0-0.14 X 5 friction: 2 nd third of process time 0-0.14 X 6 friction: 3 rd third of process time 0-0.14 X 7 hardening exponent 0.1-0.3 X 8 sheet layout 100-150% Table 6.1: Input variables of the thickness reduction application.

As the duration of the simulation lies between 2 and 6 hours per run, it is out of the scope of a direct analysis. Instead, a Latin hypercube design with 50 runs is applied as initial design for metamodeling. The design and the corresponding thickness reduction results can be found in the Appendix in Tab. B.7. It shows that the output values fall indeed into [0, 1] and are quite well spread over this interval. For the further analysis, the input variables are scaled to lie between 0 and 1. While this does not affect the sensitivity analysis, it simplifies computations and is necessary for the comparisons in the support analysis. On the results, a Kriging metamodel (2.1) with a standard product kernel, a single constant trend, and Matérn 5/2 covariance functions is fit using the R package DiceKriging [START_REF] Roustant | DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by Kriging-based metamodeling and optimization[END_REF]. Its prediction function is used as a metamodel.

The first-order Sobol and total sensitivity indices are calculated by the frequency-based methods FAST and EFAST (Section 2.3.3) using a number of n FAST = 2 000 sample points. For the computation and visualization the R package sensitivity [START_REF] Pujol | sensitivity: Sensitivity Analysis[END_REF] is used. The results are shown in Fig. 6.3. According to this, variable X 8 , the sheet layout, has the strongest influence followed by the two last parts of the friction and the flow stress. The variables X 2 , X 4 , and X 7 -the initial sheet thickness, the first third of the friction, and the hardening exponent -show only negligible influence and can thus be removed from the analysis. Nevertheless, X 4 , the first third of friction, is kept in order to be consistent in keeping friction as a parameter. The total sensitivity indices show strong interactions for all active parameters. A new Kriging model on the same data but without X 2 and X 7 is set up as updated metamodel.

X 1 X 2 X 3 X 4 X 5 X 6 X 7 X 8 0
As we are working with the very fast-to-evaluate Kriging metamodel, we can explore this main effect analysis further by looking at the support indices presented in Chapter 5. This analysis can in particular check if the chosen input ranges influence the results critically. Support indices are estimated for the six remaining variables of interest at 40 equally distributed points t = 0, . . . , 1 over the normalized domain via threshold cut 0.000 0.010 0.030 0.050 0.100 1.000 RMSE 0.793 0.585 0.780 0.735 0.786 0.860 Table 6.2: Thickness reduction application, threshold identification with standard product kernel as initial kernel, leave-one-out RMSE values for different thresholds.

In terms of interaction, the analysis so far has shown that the variables do interact, but not in which way. Thus, the specific structure of the interactions shall be examined closer using the methods presented in Chapter 3 via R package fanovaGraph. A FANOVA graph is set up with the TII estimated by the Liu and Owen method. A total number of 20 000 evaluations is used. The resulting FANOVA graph can be seen in Fig. 6.5 (b). All interactions are active, although the majority to a very small, possibly negligible, amount. The strongest interaction can be seen between X 5 and X 8 , and generally, X 5 , X 6 and X 8 seem to interact the most.

To apply the graph in subsequent analysis, a threshold to pick out the relevant interactions has to be determined. The delta jump plot (see Section 3.6) is shown in Fig. 6.5 (a). It reveals four main jumps around the values 0.01, 0.03, 0.05, and 0.1. The threshold decision procedure by Kriging kernel comparison is performed for those four threshold candidates and 0 and 1, first by using the standard product kernel as initial kernel. The resulting leave-one-out RMSE values can be seen in Tab. 6.2. The table

shows that the smallest RMSE can be achieved when removing all interactions with a TII lower than 0.01. Figure 6.5 (c) shows the corresponding FANOVA graph. The variables X 3 , the blankholder force, and X 4 , the first third of friction, are separated, while the other variables fully interact, except for X 1 and X 5 .

The threshold determination procedure is repeated with the ANOVA kernel by [START_REF] Durrande | ANOVA kernels and RKHS of zero mean functions for model-based sensitivity analysis[END_REF] as initial kernel, which showed superiority in the simulation study in Section 3.6. The results can be seen in Fig. 6 threshold cut 0.000000 0.000003 0.000200 1.000000 RMSE 0.792732 0.725182 0.723123 0.859808 Table 6.3: Thickness reduction application, threshold identification with ANOVA kernel as initial kernel, leave-one-out RMSE values for different thresholds.

delta jump plot (a) X 1 X 3 X 4 X 5 X 6 X 8 (b) X 1 X 3 X 4 X 5 X 6 X 8 (c)
tion between X 5 and X 8 being by far the strongest, indicated by a huge jump in the delta jump plot. The threshold identification procedure on Kriging kernel comparison (Tab. 6.3) even leads to an emptier graph which contains the interaction between X 5

and X 8 only.

The results of the two kernels differ clearly but at least are consistent in identifing the interaction between X 5 and X 8 as strongest. The ANOVA kernel showed superiorty in the simulation study, but in this application, the graph based on the standard product kernel lead to a smaller leave-one-out RMSE value (Tab. 6.2). To obtain more reliable results on the block-additive structure, a larger study on a new test data set seems to be necessary.

Conclusion

A complete sensitivity analysis of the thickness reduction was performed where the newly developed methods of interaction screening, support analysis, and threshold determination were applied. A strong total interaction was found between the sheet layout and the friction of the second third of the process time. Smaller interactions were identified between the sheet layout and the friction of the last third of the process time and also between both friction variables. For the identification of a thresholded FANOVA graph for subsequent studies, a deeper investigations is necessary. X 1

X 3 X 4 X 5 X 6 X 8 (b) X 1 X 3 X 4 X 5 X 6 X 8 (c)

Functional input application

Adjustments to the deep drawing press and to the simulation tools at the IUL allow the variation of the process parameters blankholder force and friction during the forming process, a possibility to further improve the shape accuracy of the formed part by more specialized parameter settings. The effect of functional friction has already been roughly examined in the thickness reduction application where it was separated in three parts considered as independent variables. This analysis shall be extended to a detailed functional sensitivity analysis of friction and blankholder force that gives insight into the functional behavior of both parameters and screens out negligible time intervals for further analysis and optimization. A small previous study showing the usefulness of the functional examination has already been shown in the motivational example in the beginning of Chapter 4. A part of the results is published in ul [START_REF] Hassan | Finite element simulations for sheet metal forming process with functional input for the minimization of springback[END_REF].

We apply the methods presented in Chapter 4 for functional sensitivity analysis to a situation with two functional inputs (friction and blankholder force) and no scalar input. This time, the scalar output is represented by the mean springback as a measure of the geometrical accuracy. Three steps of functional sensitivity analysis are performed following the design based on sequential bifurcation (Section 4.3). Mirror runs are added to ensure unbiased estimates, which results in 24 runs in total. See the Appendix for the complete design (Tab. B.8). Figure 6.7 shows the normalized regression indices of the three steps. At first, two intervals are considered for both input variables, which reveals a stronger influence of friction than blankholder force as well as a general tendency towards a positive influence in the first half and negative influence in the second half. At this point, all intervals are interesting to the engineers and thus are split for the second step. There, the last intervals of each half show to have the greatest impact whereas the first and third intervals have relatively small influence. Therefore, in the third step of the sequential analysis, these intervals are not explored further, but each second and fourth interval is split again. After these three steps, the functional sensitivities are explored sufficiently for the engineers.

To confirm the result, a further analysis using factorial designs is conducted, which requires a higher number of runs than sequential bifurcation but allows for the estimation of interactions. For the first step, a full factorial 2 4 -design with 16 runs is chosen. The second and the third step are both performed on a 2 8-3 -design with 32 runs, in which the confounding structure is chosen such that the counfounding of influential time intervals is avoided. The plot of normalized regression indices in Fig. 6.8

shows only a few dissimilarities compared to the sequential bifurcation analysis which overall confirms the chosen model. Table B.9 in the Appendix shows the estimated interaction indices together with the confounding structure. Generally, interactions are stronger for friction than for blankholder force and generally close intervals interact stronger than remote ones. The largest interaction can be observed between the intervals [11.25, 13.125] and [13.125, 15] of friction, which are directly consecutive and which are the two intervals with the strongest main effects. The second largest interaction influence is shared by the interactions between the intervals [3. 75, 5.625] and [5.625, 7.5] and the intervals [11.25, 13.125] and [13.125, 15], both corresponding to blankholder force. These intervals are again consecutive and comparatively strong.

When interpreting the general results, at first an overall similar behavior of both variables can be clearly seen, which confirms a presumption of the engineers as both parameters affect the springback in the similar way. A clear break in about half of the time can be observed. It can also be noted that the magnitude of the influence increases towards the end of both halves. These intervals are most important for the springback behavior.

Conclusion

For the essential task of springback compensation in sheet metal forming a functional sensitivity analysis was performed. To the best of our knowledge, this is the first time that the influence of time intervals is systematically investigated in sheet metal forming. From the engineering point of view, the main conclusion is that a high value of blankholder force and friction coefficient for the last 15 percent of the punch travel can strongly reduce springback in the final part. This effect can be explained by the increase in the plastic strains in the part at the end of the process, which reduces the springback. A further insight is that the parameters should have a small value during 50-60 percent of the punch travel as the sheet is being formed over the larger radius of the punch at that time. The input setting that resulted in the lowest springback of 0.47 mm, run 21 of the factorial design, is pictured in Fig. 6.9. The best run with constant settings, run 2 of the sequential bifurcation design in Tab. B.8, resulted in a springback of 0.74 mm. This reduction of the springback by around 35 percent underlines the benefits of the functional variation approach and gives an impression of its possibilities for subsequent springback optimization.

Chapter 7

Conclusion and Outlook

The general field of this thesis was sensitivity analysis for black-box functions. Three different concepts that extend the basic sensitivity analysis have been presented. They

  . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2 The role of metamodels . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3 Variance-based indices . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.3.1 Definition of various variance-based indices . . . . . . . . . . . . 2.3.2 Monte Carlo estimation . . . . . . . . . . . . . . . . . . . . . . 2.3.3 Frequency-based estimation . . . . . . . . . . . . . . . . . . . . 2.4 Derivative-based methods . . . . . . . . . . . . . . . . . . . . . . . . . 2.4.1 Morris screening by elementary effects . . . . . . . . . . . . . . 2.4.2 DGSM . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Ongoing developments in sensitivity analysis . . . . . . . . . . . . . . . 3 Sensitivity Analysis for Interaction Screening 3.1 Total interaction indices . . . . . . . . . . . . . . . . . . . . . . . . . . 3.2 TII estimation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.3 Theoretical properties of TII estimators . . . . . . . . . . . . . . . . . . 3.4 Estimating the full set of TIIs . . . . . . . . . . . . . . . . . . . . . . . 3.5 Simulation study . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . i Chapter 1

(

  2013) discover asymptotic properties of two estimators for the closed sensitivity index (2.11) scaled by the overall variance, D C I D . Both estimators use formula (2.16) for the estimation, but f 0 and the overall variance D are estimated differently. The first estimator, S 1

Figure 2 . 1 :

 21 Figure 2.1: Plot of f 1 and f 2 of Equation (2.33).

  Figure 3.1: Sensitivity analysis of the motivational example. First-order Sobol and total sensitivity indices (left), FANOVA graph of total interaction indices (right).

  or (2.21), which do not need an estimate of f 2 0 , or the estimator τ 2 proposed in Section 7 ofOwen (2013b).

  where |.| stands for the number of elements in a set and A -B = {x|x ∈ A and x / ∈ B} denotes the difference of two sets A and B.

  , and each estimate requires a number of 2(M d + L) runs. The fixing method TII estimator requires for each of the d 2 indices a number of n M C FAST estimates.

  See Tab. 3.3 for an overview of the functions used in the simulation. The functions are listed along with an identification and the distribution of their inputs. Function function 1 is taken from Fruth et al. (2014b), Ishigami from Ishigami and Homma

  3.2. Boxplots of the respective 100 estimates give an impression of the performance in terms of bias and variance for the different classes of functions and are presented in the figures 3.2 -3.4. The detailed numerical results can be found in the Appendix in Tab. B.1 -B.6.

Figure 3 .

 3 Figure3.2 shows the boxplots of the simulation results of the two test functions with low dimension and only second-order interactions. As expected in Section 3.3, negative results are observed for the RBD-FAST, Jansen, and pick-freeze estimates, but not for the Liu and Owen and the fixing method. Overall, the Liu and Owen and the fixing method TII estimator outperform the other estimators clearly, especially for influences close to zero. That means that the Liu and Owen TII estimator enables a precise detection of inactive interactions, an important task for interaction screening. The RBD-FAST TII estimator shows a rather large variance for the case with low complex functional form function 1 and a severe bias for the more complex Ishigami function.

Figure 3

 3 Figure 3.2: Simulation study on TII estimators, boxplots of 100 estimations for functions with low interactions, low number of input variables.

Figure 3 Figure 3

 33 Figure 3.3: Simulation study on TII estimators, boxplots of 100 estimations for functions with high interactions, low number of input variables.

  Figure 3.6: Delta jump plot of an artificial example (left) and corresponding FANOVA graph (right).

Figure 3 .

 3 Figure 3.7 shows a diagram of the structural setup of the package. The TII estimation is done via the function estimateGraph. If the computer experiment can be incorporated into the R environment, it can be used directly. Otherwise a metamodel can be inserted as described in Section 2.2. The function allows for the four TII estimation methods RBD-FAST, pick-freeze, Liu and Owen, and fixing method TII estimator (compare Tab. 3.1). Its input and output behavior matches with the corresponding functions in the comprehensive package for sensitivity analysis, sensitivity (Pujol et al., 2014). The package can thus easily be included into a broader study. The function estimateGraph returns an S3 object named graphlist, which contains the TII values and further information necessary to set up the FANOVA graph. Applying the generic functions print or plot to the object results in a printed output of the estimated TII values or a plot of the FANOVA graph via package igraph (Csárdi[START_REF] Csárdi | The igraph software package for complex network research[END_REF], respectively. Corresponding to Section 3.6, a threshold cut can be performed by function threshold that takes the graphlist object as input and returns a new graphlist object, whose TII values below a given threshold are set to zero and thus not printed in the FANOVA graph plot. The two methods that support the threshold decision are implemented via the functions plotDeltaJumps for the delta jump plot and thresholdIdentification for the Kriging kernel comparison. Finally, the graphlist object can be used to adapt Kriging models as described in[START_REF] Mühlenstädt | Data-driven Kriging models based on FANOVA-decomposition[END_REF] and mentioned in the motivational example. The function kmAdditive takes a data set coming from the computer experiment under study as well as the given block-additive structure in the graphlist object and returns a block-additive Kriging model by use of the popular Kriging package DiceKriging[START_REF] Roustant | DiceKriging, DiceOptim: Two R packages for the analysis of computer experiments by Kriging-based metamodeling and optimization[END_REF]. This block-additive Kriging model can then again be used to predict new data via the function predictAdditive.

  is the best constant: If equalities can be achieved in the Poincaré inequalities satisfied by each distribution, then the inequality (3.14) reaches equality as well.Two different proofs are given in[START_REF] Roustant | Crossed-derivative based sensitivity measures for interaction screening[END_REF], one via sequential application of Poincaré inequalities and one using the crossed finite difference in the Liu and Owen formula (3.5).

Figure 4 . 1 :

 41 Figure 4.1: Preview runs. Comparison of different functional friction settings.

Figure 4 . 2 :

 42 Figure 4.2: Visualization of the piecewise constant representation of functional inputs.

  be estimated using a design matrix Z on the values of Z (k) j sufficient for estimation, i.e. such that Z Z is invertible (see e.g.[START_REF] Saltelli | Wiley series in probability and statistics[END_REF], p. 124). The estimates α and β (k) j are obtained by the Least Squares formula (Z Z) -1 Z y, with y the vector of the output realizations. As the underlying computer experiment is deterministic, no assumptions are made about ε, the difference between response and regression model and the Least Squares approach is interpreted as simple curve fitting. The coefficients β (k) j can then be used as sensitivity indices as they indicate the linear influence of the corresponding functional interval on the output.

  Figure 4.3 (top right) shows this situation. In contrast to the previous plot, the plots indicates an increasing influence towards the end of the functional domain.Independence of the chosen decomposition can be achieved by normalization through the interval size.Definition 4.1. Consider a set of splitting points a 1 , . . . , a p d and assume that g j ∈ V a j , j = 1, . . . , d : g j (t) = Z

Figure 4 . 3 :

 43 Figure 4.3: Graphical representation of the functional sensitivity of an artificial example, via regression coefficients β (top and bottom left, for different splittings) and via normalized regression coefficients H (bottom right). The bar colors highlight size and direction of the indices. Unnormalized, the regression coefficients are sensitive to the splitting.

Figure 4 . 4 :

 44 Figure 4.4: Scheme for the design based on sequential bifurcation. White and gray shading indicate the setting to +1 and -1.

  Fig. 4.5 for an overview of the package structure. Its core is an object, safidesign. It contains the current design, which can be enhanced to get to the next step of the sequential procedure. The splitting can be performed following the sequential bifurcation algorithm or by any desired design given manually. Mirror runs can be added. At any step, the object can be accessed to obtain a transferable design matrix, and also be plotted. A model can be fit to corresponding output from the computer experiment. The resulting safimodel object contains the normalized regression indices computed according to the chosen design. Applying the function plot to this object leads to the graphical representation by barplots as shown in Fig. 4.3.

  Figure 4.5: Structural setup of R package seqSAFI. Functions are represented by rectangles (with double vertical lines if defined externally), data by parallelograms, and output by circles. Arrows represent information transfer.

  Fig. 4.6. The six runs are executed via Equation (4.4) and the normalized regression indices of Def. 4.1 are computed as described in Section 4.3. They are visualized in the bottom of Fig. 4.6.The first functional input g 1 shows only negligible influence in its second interval. Thus, only the first interval shall be explored further and is split into three intervals. For g 2 , no intervals can be screened out, so that both intervals are split up. This results in four new runs depicted in the top of Fig.4.7. The resulting indices can be seen in the bottom of the same figure.

Figure 4 . 6 :

 46 Figure 4.6: Analytical example, step 1, design and normalized regression indices.

Figure 4 . 7 :

 47 Figure 4.7: Analytical example, step 2, design and normalized regression indices.

Figure 4 . 8 :

 48 Figure 4.8: Analytical example, step 3, design and normalized regression indices.

=

  Minimum and maximum of the input variable distribution are usually specified by the practitioner as extreme values of the parameter in question. Those values are often chosen vaguely, but can have an important impact. As an example, the total effects of the Ishigami function with the usual distributionsX i ∼ U [-π, π] are D T 1 = 7.72, D T 2 = 6.13, and D T 3 = 3.37. If we reduce the support of the distribution of X 3 by 10 percent, X * 3 ∼ U [-π + π /10, π -π /10], we get D T * 1

  Consider again the standard situation of a black-box function Y = f (X), where X = (X 1 , . . . , X d ) is a vector of independent random variables with distribution µ= µ 1 ⊗ • • • ⊗ µ d and f a d-dimensional function f : ∆ → R. In addition we assume that ∆ = [0, 1] d compact, f of class C 1 ,and that X 1 , . . . , X d have continuous density functions with support [0, 1]. Notice that since f is of class C 1 on the compact set [0, 1] d , f and all its first-order derivatives are bounded, ensuring that f (X) and ∂f ∂x• (X) are in L 2 .

  Now remark that R is a bounded random variable by continuity of ∂ 2 f ∂x 2 i on the compact set [0, 1] d , and thus R is bounded as well. This implies that Var(R) = O(1) and Cov(B i , R) = O(1), and the result follows.

  Figure 5.1: Support analysis of the Ishigami function.

FigFigure 6

 6 Fig. 6.1.

  Figure 6.3: Thickness reduction application, first-order Sobol and total sensitivity indices.

Figure 6 . 5 :

 65 Figure 6.5: Thickness reduction application with standard product kernel as initial kernel. Delta jump plot (a). FANOVA graphs of the complete interaction structure (b) and the thresholded values (c).

Figure 6 . 6 :

 66 Figure 6.6: Thickness reduction application with ANOVA kernel as initial kernel. Delta jump plot (a). FANOVA graphs of the complete interaction structure (b) and the thresholded values (c).

Figure 6 . 7 :

 67 Figure 6.7: Functional input application, normalized regression indices of three steps (top, middle, bottom). The bar color references sign and amplitude of the bar.

Figure 6 . 8 :

 68 Figure 6.8: Functional input application in three steps via factorial designs. First step (top). Second step (middle). Third step (bottom). The bar color references sign and amplitude of the bar.

  all share the ability to give deeper insight into the function's behavior on different levels, be it interactions, the domain of functional inputs, or the distribution support of inputs.The sensitivity analysis for interaction screening enables a user to explore the interaction structure of the function and to find inactive interactions. This provides a block-additive decomposition of the black-box function, which in turn can be of use, for instance in prediction modeling and optimization.The sensitivity analysis for functional inputs gives insight into the influence of the inputs on the functional domain level. The method explores the influence of different parts of the functional domain, possibly leading to enhanced insights about the influence of the input variables compared to scalar analysis. The screening of inactive parts of the functional domain can again serve as a starting point for modeling and optimization.The support analysis, finally, provides an extension from single sensitivity indices to sensitivity functions by returning a sensitivity for each point of the distribution support of the input. An insight into the influence of the inputs for different supports can be gained, showing the potential consequences of input distribution choices. symbol description mathematical notations L p (µ) space of functions that are p times integrable with respect to measure µ C k class of functions with continuous derivatives up to order k x Euclidean norm over vector x ∇f (X) gradient of f |Q| for a set Q, number of elements in Q A -B {x|x ∈ A and x / ∈ B}, difference of sets A and B F [0,1] space of functions with domain [0, 1] C(µ) Poincaré constant of a measure µ defined variable names (general) n size of Monte Carlo sample d number of input variables x = (x 1 , . . . , x d ) vector of input setting (for one run) X = (X 1 , . . . , X d ) random vector input variables k(h) Kriging covariance kernel for input distance h (2.3) δ * small real-valued spacing, e.g. in finite differences I set of variable indices, subset of {1, . . . , d} x I , X I input setting and variable for all input variables in I x -I , X -I input setting and variable for all input variables but the ones in I y, Y specific output value and corresponding random variable f underlying function, either computer experiment or metamodel µ = µ 1 ⊗ • • • ⊗ µ d probability measure of X ∆ domain of f f I additive term of f in FANOVA decomposition (2.4) in Chapter 4 d scal , d fun number of scalar and functional input variables g 1 , . . . , g d fun functional input variables p j number of intervals of the domain decomposition of input g j input g j over interval [a k-1 , a k ] (4

Table 3

 3

.1: Overview of considered TII estimators. where D k |X -{i,j} denotes the overall variance of f fixed and D k i|X -{i,j} and D k j|X -{i,j}

Table 3 .

 3 2: Number of function evaluations for the considered TII estimators.

Table 3 .

 3 3: Overview of test functions.

		d ≤ 5		d ≥ 10
		form = simple form = complex
	interaction = low	function 1	Ishigami	high 2 nd order
	interaction = high pure 3 rd order	Branin4	high 14 th order

Table 3 .

 3

4: Classification of test functions.

Acknowledgements

finite differences with a Monte Carlo sample of size 2 000 for each point. The resulting support functions are drawn in Fig. 6.4. The plots show that for the most influential variables X 8 , X 6 , X 5 , and X 1 the specific range indeed does not affect the result strongly, as the output variation is not strong at the borders of the input support. The wavy shape of the support index function of variable X 8 reveals strong nonlinearities in the relation between this variable and the thickness reduction. The two friction parameters X 5 and X 6 show similar support functions, underlining that both describe the same parameter. For X 3 and X 2 , the function of the input influence over the domain does not change much, possibly due to their small overall influence. For all variables, the interactions seem to be rather equally spread over the support.

Conclusion and Outlook

The presented methods were applied in two forming process studies, concerning the thickness of the formed metal and the springback, respectively. They helped to improve the understanding of the forming process and gave more insight into the way input sources influence the deformation yielding an important prerequisite to subsequent research in the field of sheet metal forming.

All three methods provide complete procedures for the exploration of computer experiments, which can lead to new insights and can substantially support further analyses. There are, nevertheless, several possibilities for improvement. For the sensitivity analysis of interaction screening, several extensions are imaginable, e.g. to allow for the employment of dependent variables or to explore metamodels that enable the direct computation of the TIIs, possibly via ANOVA zero-mean kernels or Polynomial Chaos Expansion. In the sensitivity analysis for functional input, aims for future research could be to develop methods that allow for smooth design functions, non-linear influences, and models that avoid the problem of canceling out variables of opposite sign. However, those will be challenging tasks due to the necessary strong dimension reduction and the focus on sensitivity analysis. One further point of improvement is to consider not only the value of the input at the time points but also how the values between the time points evolve, and thus to set up models that incorporate this evolution in time. A straightforward outlook is the extension to spatial input variables.

One idea could be to treat each spatial dimension as a functional input and thus lead the problem back to functional analysis. Another idea is to split the input space into simplices similar to the splitting of intervals in the functional analysis.

Appendix