Prof Olivier Le Maître

Prof. Jérôme Morio

Prof. Benoît Robyns

Prof Prof Marina Kleptsyna

Michel

Anna Kharlamova

Fan Wang

Lei Xiong

Liyang Zheng

Lizhong Zhao

Ru Wang

Van Tang Nguyen

Keywords: Uncertainty analysis, Sensitivity analysis, Variance decomposition, Sampling, Correlation, Sensitivity measure, Opinion dynamics i

By means of taylor series expansion, a general analytic formula is derived to characterise the uncertainty propagation from input variables to the model response, in assuming input independence. By using power-law and exponential functions, it is shown that the widely used approximation considering only the first order contribution of input uncertainty is sufficiently good only when the input uncertainty is negligible or the underlying model is almost linear. The method is then applied to a power grid system and the EOQ model.

The method is also extended to correlated case. With the extended method, it is straightforward to identify the importance of input correlations in the model response.

This allows one to determine whether or not the input correlations should be considered in practical applications. Numerical examples suggest the effectiveness and validation of our method for general models, as well as specific ones such as the deterministic HIV model.

Our method is then compared to Sobol's one which is implemented with sampling based strategy. Results show that, compared to our method, it may overvalue the roles of individual input factors but underestimate those of their interaction effects when there are nonlinear coupling terms of input factors. A modification is then introduced, helping understand the difference between our method and Sobol's one.

Finally, a numerical model is designed based on a virtual gambling mechanism, regarding the formation of opinion dynamics. Theoretical analysis is proposed by the use of one-ata-time method. Sampling-based method provides a global analysis of output uncertainty and sensitivity.

Résumé

Par un développement en série de Taylor, une relation analytique générale est établie pour calculer la propagation des incertitudes des variables d'entrée sur la réponse du modèle, en assumant l'indépendance des entrées. En utilisant des relations puissances et des relations exponentielles, il est démontré que l'approximation souvent utilisée consistant à ne considérer que la contribution du premier ordre sur l'incertitude d'entrée permet d'évaluer de manière satisfaisante l'incertitude sur la réponse du modèle pourvu que l'incertitude d'entrée soit négligeable ou que le modèle soit presque linéaire. La méthode est appliquée à l'étude d'un réseau de distribution électrique et à un modèle d'ordre économique.

La méthode est étendue aux cas où les variables d'entrée sont corrélées. Avec la méthode généralisée, on peux déterminer si les corrélations d'entrée doivent ou non être considérées pour des applications pratiques. Des exemples numériques montrent l'efficacité et la validation de notre méthode dans l'analyse des modèles tant généraux que spécifiques tels que le modéle déterministe du VIH.

La méthode est ensuite comparée à celle de sobol. Les résultats montrent que la méthode de sobol peut surévaluer l'incidence des divers facteurs, mais sous-estimer ceux de leurs interactions dans le cas d'interactions non linéaires entre les paramètres d'entrée. Une modification est alors introduite, aidant à comprendre la différence entre notre méthode et celle de sobol.

Enfin, un modéle numérique est établi dans le cas d'un jeu virtuel prenant en compte la formation de la dynamique de l'opinion publique. L'analyse théorique à l'aide de la méthode de modification d'un paramètre à la fois. La méthode basée sur l'échantillonnage fournit une analyse globale de l'incertitude et de la sensibilité des observations. Mots-clés: Analyse d'incertitude, Analyse de sensibilité, Décomposition de variance, Echantillonnage, Corrélation, Mesure de sensibilité, Dynamique d'opinion First and foremost, I would like to show the greatest appreciation to my supervisor, Prof. Q. Alexandre Wang, who provided me an exciting opportunity to study in France, particularly, in ISMANS at the beginning and then at Le Mans University. He guided me to a more international research environment and also offered me a special platform to communicate and discuss with scholars from different fields and countries. It has been a great honor to be one of his Ph.D. students. I deeply appreciate his continuous support, innovative ideas, valuable suggestions and timely guidance, which all are useful and beneficial to my Ph.D. study and related research, and also make my study experience productive and stimulating. The tremendous interests and enthusiasm he has for his research, as well as his immense knowledge, are contagious and motivational for me, even during tough times in the Ph.D. pursuit. I am also grateful for the great help that he and his wife offered for my daily life in France.

I am especially thankful to my supervisor, Prof. Xu Cai, working at Central China Normal University. He introduced me to the research filed of Complexity Science and provided me many chances to attend international academic conferences. He encouraged me to persistently move forward and to work independently. As a mentor, except for the knowledge related to research, he also taught me how to positively face setbacks in study and life, how to communicate with others, how to build good personality, how to "conduct oneself", etc. His encouragement and unlimited support incented me to do better both in study and life.

My sincere thanks also go to Prof. Wei Li, my another supervisor working at Central China Normal University. He guided me to do numerical simulations and theoretical analysis. He taught me how to write a research paper, how to reply the referee's report, etc. I am particularly indebted to him for his constant encouragement and endless support throughout the whole stage of my postgraduate study. His insightful comments and valuable ideas have particularly contributed to the writing of this thesis. I would like to thank my associate supervisor, Prof. Alain Bulou, for his constant encouragement and support when I studied at IMMM. He provided me financial support to attend a particular international academic conference, and permission to use facilities of the laboratory. He is very gentle and always offers help to me when I am in troubles.

iii

General introduction

The concept of global uncertainty and sensitivity analysis has proposed for a long time.

It attracts the considerations of a large number of researchers from various of areas. The global uncertainty and sensitivity analysis aims at analysing the uncertainties of output variables (also called observations or model responses) according to the uncertainties in input variables (or named factors, parameters, covariates), and the sensitivity of each output variable with respect to individual input parameters, as well as to their interactions. Undoubtedly, global uncertainty and sensitivity analysis is advantageous for gaining insight into how input variables should be ranked according to their importance in establishing the uncertainties in different output variables. While the strategies for uncertainty and sensitivity analysis are quite extensive, a general analytic method is still limited, especially for models of present input correlations. In this dissertation, we mainly focus on the establishment of a general theoretical framework for global uncertainty and sensitivity analysis in the modelling of complex systems.

Mathematical models have a wide range of applications in diverse disciplines. They can help explain a system by visualized data and/or figures and analyse the possible effects of different parameters, and also, if necessary, make predictions about the underlying behaviour. With a deterministic mathematical model of general form y = f (x) with y denoting the output vector and x indicating the input vector. When y is calculated from x through a specified function or some natural or artificial rules connecting y and x, uncertainties in the elements of input vector, if exist, will propagate through the calculation to the members of output vector y dependent on x [1,2]. This process is called variance propagation (or uncertainty propagation). Variance propagation, often regarded as the fundamental ingredient of sensitivity analysis for complex models, mainly considers the determination of output's variance via uncertainties in input factors [3,4].

At present, many strategies have been built for the determination of variance propagation, such as simulation-based methods [5,6], most probable point-based methods [7,[START_REF] Fiessler | Quadratic limit states in structural reliability[END_REF], functional expansion-based methods [START_REF] Xiu | Modeling uncertainty in flow simulations via generalized polynomial chaos[END_REF], numerical integration-based methods [START_REF] Seo | Efficient statistical tolerance analysis for general distributions using three-point information[END_REF][START_REF] Xu | A generalized dimension-reduction method for multidimensional integration in stochastic mechanics[END_REF][START_REF] Rahman | A univariate dimension-reduction method for multidimensional integration in stochastic mechanics[END_REF][START_REF] Lee | Response surface augmented moment method for efficient reliability analysis[END_REF].

Simulation-based methods, also called sampling-based methods, are regarded as both effective and widely used, especially for those models with the functional relationship connecting y and x absent [START_REF] Blower | A tale of two futures: Hiv and antiretroviral therapy in san francisco[END_REF][START_REF] Heltona | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF][START_REF] Marcot | Analysis of sensitivity and uncertainty in an individual-based model of a threatened wildlife species[END_REF]. These briefly mentioned strategies, however, are computationally expensive, especially in the presence of a high number of input variables. For a general model with given functional form, the procedure will be much easier and numerically cheaper for determining the output's variance if an analytic formula associated with variance propagation can be provided. More information associated with other methods for variance propagation can be found in some reviewed papers [START_REF] Padulo | Comparative analysis of uncertainty propagation methods for robust engineering design[END_REF][START_REF] Lee | A comparative study of uncertainty propagation methods for black-box-type problems[END_REF][START_REF] Borgonovo | Sensitivity analysis: A review of recent advances[END_REF].

A simple analytic formula has been appeared since 1953. It approximately computes the variance of the product of two independent random variables [START_REF] Yates | Sampling Methods for Censuses and Surveys[END_REF]. In 1966, this approximation was extended by engineers and experimentalists to more general multivariate cases [START_REF] Ku | Notes on the use of propagation of error formulas[END_REF]. This formula, also called Taylor series approximation, restricted to the first-order terms [START_REF] Helton | Uncertainty and sensitivity analysis techniques for use in performance assessment for radioactive waste disposal[END_REF], has gained a wide applications thanks to its simplicity and convenience [START_REF] White | Uncertainty of measurement in quantitative medical testing a laboratory implementation guide[END_REF]. However, it can satisfactorily estimate the output's variance only when the functional relationship between output and input variables is almost linear or the uncertainty of each input variable is negligible [START_REF] Padulo | Comparative analysis of uncertainty propagation methods for robust engineering design[END_REF]. For most models, y highly nonlinearly depends on x having large uncertainties. This suggests the necessity of an exact analytic formula in calculating the output's variance and evaluating its sensitivities with respect to individual input factors, as well as to their interactions.

Furthermore, many methods have also been designed for performing sensitivity analysis, including the traditional approach of changing one factor at a time [START_REF] Guo | Heterogeneity and subjectivity in binary-state opinion formation systems[END_REF][START_REF] Zhu | Opinion evolution on a ba scaling network[END_REF], local method [START_REF] Cacuci | Sensitivity and uncertainty analysis: Applications to Large-Scale Systems[END_REF], regression analysis [START_REF] Iooss | A review on global sensitivity analysis methods[END_REF], variance-based method [START_REF] Zhu | An analytic method for sensitivity analysis of complex systems[END_REF], etc. Among the various available strategies, variance-based sensitivity analysis has been assessed as versatile and effective for uncertainty and sensitivity analysis of model response. The consideration of variance-based importance measures can be traced back to over twenty years ago when Sobol characterised the first-order sensitivity measures on the basis of deposing the variance in model response into different partial contributions attributable to individual input variables and to their combinations (called variance decomposition) [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. Then extensive relevant investigations are carried out around this Sobol's work, boiling down to the improvements in analysis strategies and to their applications to the sensitivity and reliability analysis of complex systems [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF][START_REF] Kucherenko | Comparison of lagin hypercube and quasi monte carlo sampling techniques[END_REF]. However, these frameworks, as well as above mentioned strategies for the determination of variance propagation, are often proposed when the input variables are assumed to be statistically independent.

Recently, the interest in extending sensitivity analysis strategies from uncorrelated case to the correlated one is increasing as correlated factors are often happened in practical applications. Previous investigations about sensitivity analysis of models in the presence of input correlations only provided overall sensitivity indices with respect to individual input factors. However, the correlated and independent variance contributions were Chapter 1

Introduction

Mathematical models are of great importance in the natural sciences. They have been diffusely utilized in many disciplines as diverse as mathematics [START_REF] Kermack | A contribution to the mathematical theory of epidemics[END_REF][START_REF] Mardia | Multivariate analysis (probability and mathematical statistics[END_REF], physics [START_REF] Guo | Heterogeneity and subjectivity in binary-state opinion formation systems[END_REF][START_REF] Zhu | Opinion evolution on a ba scaling network[END_REF][START_REF] Castellano | Statistical physics of social dynamics[END_REF][START_REF] Kitsak | Identification of influential spreaders in complex networks[END_REF][START_REF] Pei | Spreading dynamics in complex networks[END_REF][START_REF] Hui | Structure properties of evolutionary spatially embedded networks[END_REF], chemistry [START_REF] Sterratt | Principles of computational modelling in neuroscience[END_REF], etc. With mathematical models, one can explicate a system in mathematical language, analyse the roles of linked factors by physical methods, and then make reasonable predictions of underlying behaviors. In general a model contains three major elements: the input vector, the output vector, and associations between them. In practical applications, the elements of input vectors are rarely deterministic but contain uncertainty following some distribution laws [START_REF] Vaughn | Representation of two-phase flow in the vicinity of the repository in the 1996 performance assessment for the waste isolation pilot plant[END_REF][START_REF] Campolongo | Hitchhiker's guide to sensitivity analysis[END_REF]. Consequently, the determination of the variations in input variables, the investigation of their propagating through the model, as well as the quantification of the sensitivities of model outputs with respect to input variables are of crucial importance for establishing reliable and robust models [3,[START_REF] Xu | Uncertainty and sensitivity analysis for models with correlated parameters[END_REF][START_REF] Helton | Uncertainty and sensitivity analysis for models of complex systems[END_REF][START_REF] Saltelli | Global sensitivity analysis: the primer[END_REF]. The implementation of these procedures is known as uncertainty and sensitivity analysis.

A view of modeling that may help illustrate the role of uncertainty and sensitivity analysis in the scientific process is offered in Fig. 1.1, taken from the work of Robert Rosen, an American theoretical biologist [START_REF] Rosen | Life itself: a comprehensive inquiry into the nature, origin, and fabrication of life[END_REF]. The figure shows two systems, a natural system N which forms the subject under investigation, and a formal system F which indicates the modeling of this subsect. Each system has its own internal entailment structures and the two systems are connected by the encoding and decoding processes.

The uncertainty under discussion here is often referred to as epistemic uncertainty (also known as systematic uncertainty). Epistemic uncertainty derives from a lack of information or non-accuracy in measurement about the appropriate value used for specifying a quantity that is assumed to be constant in the context of the analysis for a particular problem. In the conceptual and computational designation of an analysis, epistemic uncertainty is regarded in general to be distinct from aleatory uncertainty, which, also known as statistical uncertainty, arises from an inherent randomness in the behavior of the system under study.

Uncertainty and sensitivity analysis are essential parts of analyses for complex systems.

Specifically, sensitivity analysis considers the determination of variance contributions of individual input variables to the elements of output vectors [START_REF] Saltelli | Sensitivity analysis in practice: a guide to assessing scientific models[END_REF]. Uncertainty analysis, preceding sensitivity analysis, mainly focuses on the determination of the uncertainties in output variables that derive from uncertainties in input factors. Conceptually, uncertainty and sensitivity analyses should be run in tandem. They work together to help determine:

(1) Which input factors contribute most to the variation of model output.

(2) Which parameters are significant and which ones can be eliminated from the model;

(3) How to efficiently reduce the uncertainty in model output by strengthening the knowledge base concerning input parameters.

Quantifying the impact of a variable under sensitivity analysis could be useful for a series of purposes, such as deep understanding of the relationships between input and output variables in a system or model, fixing model inputs that have no effect on the output, identifying and removing redundant parts of the model structure, and avoiding useless time consumption on non-sensitive variables in models of a large number of parameters. In models consisting of a large number of input variables, sensitivity analysis constitutes essential ingredient of model building and quality assurance. Sensitivity analysis has also extended its application to national and international agencies involved in impact assessment studies, including the European Commission [50, 51], Australian pathology laboratories [START_REF] Farrance | Uncertainty of measurement: a review of the rules for calculating uncertainty components through functional relationships[END_REF], the Intergovernmental Panel on Climate Change [START_REF]Intergovernmental Panel on Climate Change. Climate Change 2014-Impacts, Adaptation and Vulnerability: Regional Aspects[END_REF], and US Environmental Protection Agency's modelling guidelines [START_REF] Gaber | Guidance on the development, evaluation, and application of environmental models[END_REF].

The framework of uncertainty and sensitivity analysis is easily performed when only a single input factor is involved in the model under discussion, which is known as univariate situation. It requires a straightforward one-dimensional analysis by presenting results in figures in a two-dimensional space. When two or more input factors are under assessment, however, the problem is much more complicated, especially if input factors do not have a separable monotonic effect on the output variable of interest.

The early framework of sensitivity analysis for multivariate models was established by using local analysis. Local sensitivity analysis aims at assessing the local effects of uncertainties in individual input factors on output variables, by concentrating on the sensitivity in vicinity of a set of special factor values [START_REF] Gustafson | Local sensitivity analysis[END_REF]. Such sensitivity is usually evaluated by the use of gradients or partial derivatives of functions connecting output and input vectors at these special factor values. This means the values of the rest factors are fixed while studying the local sensitivity of model response with respect to a single factor.

Local sensitivity analysis is most frequently employed for the analysis of complex models, especially when a large number of input factors are involved. This is common because of the simplicity and low computational cost in its implementation. However, it abortively quantifies the global impact of individual input factors and of their interactions on output variables. Of importance to a part of model analysis practitioners (mostly working in the fields of statistics, risk and safety assessment, and reliability detection) is understanding the sensitivity of an output variable with respect to simultaneous variations of several input factors [START_REF] Saltelli | Global sensitivity analysis: the primer[END_REF][START_REF] Wagner | Global sensitivity analysis[END_REF]. The global uncertainty and sensitivity analysis provides such sensitivity information. It evaluates the influence of individual input factors by looking at the entire input space rather than at a specified point.

Generally, the process of global uncertainty and sensitivity analysis can be decomposed into: (a) specifying the model under study and defining its input and output variables;

(b) characterising the uncertainty in input variables; (c) determining the uncertainty in model output; and (d) quantifying the importance of individual input variables in the estimation of output variables. In establishing the framework of sensitivity analysis for a given model with defined input and output variables, the main goal is to handle the remaining procedures.

A model can be stated as diagnostic or prognostic. Diagnostic models are used for understanding a law. They are often built by wild speculations applied to play what-if games, such as models designed to study the emergence of an agreement in a population [START_REF] Maity | Emergence of fast agreement in an overhearing population: The case of the naming game[END_REF][START_REF] Beuls | Agent-based models of strategies for the emergence and evolution of grammatical agreement[END_REF][START_REF] Centola | The spontaneous emergence of conventions: An experimental study of cultural evolution[END_REF], models used to investigate organizational change, etc. In the investigation of organizational change, as an example, three diagnostic models are high potential candidates to highlight the problem areas and provide structure for solution development. The first one is an analytic model, also known as the difference-integration model. It focuses on thorough analytical diagnosis as the foundation for organisational change [START_REF] Cummings | Organization development and change[END_REF]. The second one is the force-field analysis model, originally developed by Kurt Lewin in the early 1950s. It regards the organisation as the result of internal forces that drive change or maintain the current status [START_REF] Lewin | Field Theory in Social Science: Selected Theoretical Papers[END_REF]. The third one is developed on the bases of cause maps and social network analysis. It provides a mathematical approach to organisation diagnosis [START_REF] Nelson | Cause maps and social network analysis in organizational diagnosis[END_REF]. Regarding diagnostic models, input variables are pre-defined by model designers and often assumed to follow particular distribution laws in fixed real ranges (e.g. uniform and Gaussian distributions). Prognostic models can be viewed as accurate and trusted predictors of a system. They mainly focus on the estimation (prediction) of the probability that a particular event or outcome will happen. Prognostic models are often developed for the clinical practice, where the risk of disease development or disease outcome (e.g. recovery from a specific disease) can be calculated for individuals by combining information across patients. In the clinical practice, prognostic models can be presented in the form of a clinical prediction rule [START_REF] Harrell | Multivariable prognostic models: issues in developing models, evaluating assumptions and adequacy, and measuring and reducing errors[END_REF][START_REF] Altman | What do we mean by validating a prognostic model?[END_REF][START_REF] Egger | Prognosis of hiv-1-infected patients starting highly active antiretroviral therapy: a collaborative analysis of prospective studies[END_REF][START_REF] Moons | Prognosis and prognostic research: application and impact of prognostic models in clinical practice[END_REF]. Prognostic models also find their applications in other fields, such as risk and safety evaluation in engineering [START_REF] Williams | Safety evaluation and risk assessment of the herbicide roundup and its active ingredient, glyphosate, for humans[END_REF][START_REF] Klinke | A new approach to risk evaluation and management: risk-based, precaution-based, and discourse-based strategies[END_REF], problems solving of water dynamics in estuaries [START_REF] Kochergin | Three-dimensional prognostic models[END_REF]. It is often preferable that input variables in prognostic models, in contrast to those in diagnostic models, are easily determined by social experience or practical examples for ensuring the applicability of a prognostic model in practical applications. Furthermore, models can also be classified as data-driven or law-driven. A data-driven (or inverse) model tries to derive properties statistically by empirical study. Advocates of data-driven models like to describe social behaviours with a minimum of adjustable parameters, for instance, models helping understand the spreading of really happened epidemics [START_REF] Riley | Transmission dynamics of the etiological agent of sars in hong kong: impact of public health interventions[END_REF][START_REF] Lipsitch | Transmission dynamics and control of severe acute respiratory syndrome[END_REF][START_REF] Alam | A model for hiv spread in a south african village[END_REF][START_REF] Beck | A data-driven simulation of hiv spread among young men who have sex with men: the role of age and race mixing, and stis[END_REF], models designed for explaining the generation of traffic jams [START_REF] Smith | Comparison of parametric and nonparametric models for traffic flow forecasting[END_REF][START_REF] Orosz | Traffic jams: dynamics and control[END_REF], and models proposed to describe financial time series [START_REF] Hassan | Stock market forecasting using hidden markov model: a new approach[END_REF][START_REF] Dablemont | Modelling and forecasting financial time series of" tick data" by functional analysis and neural networks[END_REF]. Law-driven (or forward) models, on the other hand, aim at employing appropriate laws which have been attributed to the system to predict its behaviour. For example, people can use Darcy's and Fick's laws to understand the motion of a solute in water flowing through a porous medium [START_REF] Hassanizadeh | Derivation of basic equations of mass transport in porous media, part 2. generalized darcy's and fick's laws[END_REF][START_REF] Whitaker | The equations of motion in porous media[END_REF]. In building design, as another example, building energy simulation models are generally classified as prognostic law-driven models by which the behaviour of a complex system can be predicted in terms of a set of well-defined laws (e.g., mass balance, energy balance, conductivity, and heat transfer, etc.) [START_REF] Coakley | A review of methods to match building energy simulation models to measured data[END_REF]. In datadriven models, input parameters are introduced based on special situations, which could be deterministic and attached with an artificially defined uncertainty. For law-driven models, input factors are often imported by the laws we employed. Characterising the uncertainty in input factors is also dependent upon the situation under analysis.

With given distribution laws, the uncertainty in input variables can be specified by their mean values (e.g., arithmetic mean or mathematical expectation, geometric mean, median), standard deviations, PDFs, CDFs and CCDFs. Particularly, when input variables are deterministic (often happening in agent-based systems where input parameters can be determined based on practical experience), their uncertainties are frequently represented by artificially introducing fixed variation around their normative values or few typical scenarios (e.g., scenarios corresponding to any possible combinations of specified low, medium, high values of input factors) for performing the uncertainty and sensitivity analysis of the system under discussion [START_REF] Marcot | Analysis of sensitivity and uncertainty in an individual-based model of a threatened wildlife species[END_REF][START_REF] Happe | Agent-based modelling and sensitivity analysis by experimental design and metamodelling: an application to modelling regional structural change[END_REF][START_REF] Perez | An agent-based approach for modeling dynamics of contagious disease spread[END_REF][START_REF] Fonoberova | Global sensitivity/uncertainty analysis for agent-based models[END_REF][START_REF] Usman | Sensitivity analysis of the parameters of a model of hiv-1 transmission[END_REF]. The analysis framework of deterministic situations is often designed according to the variation in model output driven from the independent variation in each input factor. This method is known as one-at-a-time method and will be discussed below in detail. Examples, as shown in Fig. 1.2, present the characterisation of uncertainties in input factors for both kinds of situations, in the absence of input correlations.

Presentation of uncertainty in model output

As already mentioned at the beginning of this chapter, uncertainty and sensitivity analysis should be run in series, with uncertainty analysis preceding in current practice.

Uncertainty analysis is, through a certain way, to determine the uncertainty in model output based on the uncertainty in model input.

One popular way of establishing uncertainty analysis is dependent upon the computer, which is also known as Monte Carlo (MC) method. Consider a general model of the

form y = f (x) with x = (x 1 , x 2 , • • • , x n ) T indicating an input vector of n-dimensional
variables. All elements of input vector are assumed to be independent of each other.

By given PDFs of individual input factors, a sample of size M , indicated by an M × n matrix, can be generated as

          x 1 1 x 2 1 • • • x n-1 1 x n 1 x 1 2 x 2 2 • • • x n-1 2 x n 2 . . . . . . . . . . . . . . . x 1 M -1 x 2 M -1 • • • x n-1 M -1 x n M -1 x 1 M x 2 M • • • x n-1 M x n M           . (1.1)
Run independently the model for all points that are sampled in n-dimensional input space. A set of values of the model output y are then generated accordingly: The ends of the whiskers can represent several possible alternative values including: the minimum and maximum of all of the data, the 9th percentile and the 91st percentile, the 2nd percentile and the 98th percentile [START_REF] Frigge | Some implementations of the boxplot[END_REF][START_REF] Hubert | An adjusted boxplot for skewed distributions[END_REF]. In practical applications, the model output is not always a scalar but could also be a function. For example, in the investigation of epidemic spreading, the system of interest is time-dependent. Uncertainty in input factors will be propagated to the uncertainty in the dependence of model output upon the time parameter. For such situations, an effective presentation format of the uncertainty in model output is to use two graphical frames, with first one displaying any possible dependence of model output upon a relevant parameter and second one presenting statistical results for the outcomes in the first one [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF]. input factors.

y = (y 1 , y 2 , • • • , y M ) T . (1.

Methods of sensitivity analysis

Sensitivity analysis is a primary part of model development. It involves importance evaluation of input parameters in the estimation of the model output. At present, a large number of approaches have been built for performing the sensitivity analysis.

Regarding the complexity of models, many methods are developed to address one or more constrains. For example, most common sensitivity analysis methods assume input factors are independent of each other [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF][START_REF] Saltelli | A quantitative model-independent method for global sensitivity analysis of model output[END_REF]; approaches based on linear regression are valid only for linear models; virtually all sensitivity analysis methods consider a single univariate model output, by which sensitivity measures are hard to be interpreted for models with correlated outputs.

Methods of sensitivity analysis can be classified based on the methodology as mathematical, statistical or graphical [START_REF] Frey | Identification and review of sensitivity analysis methods[END_REF]. Mathematical methods evaluate the sensitivity of model output with respect to the range of variation of each input factor. Typically, they involve the calculation of output variable according to a few values of each factor that represent the possible variation range of the factor. Mathematical methods can identify the influences of individual factors in their variation ranges on an output variable. However, they do not indicate the variance of output variable propagated from the uncertainties in input factors but represent, for example, the sensitivity of model output as the magnitude of percentage change compared to its nominal value. In some cases, especially when input factors are almost deterministic, mathematical methods are helpful in recognising the most important factors [START_REF] Happe | Agent-based modelling and sensitivity analysis by experimental design and metamodelling: an application to modelling regional structural change[END_REF][START_REF] Benke | Parameter uncertainty, sensitivity analysis and prediction error in a water-balance hydrological model[END_REF]. One-at-a-time method (discussed below) is one of most widely used mathematical methods. Statistical methods assess the variance contribution of input factors to the output variable with given probability distributions related to input factors. By employing statistical methods, the variation of one or more input factors can be considered simultaneously. This allows one to identify the interaction effects among multiple input factors on the uncertainty in output variable. Some statistical methods often considered are introduced here, including differential-based method, response surface method, regression analysis, variance-based analysis, and sampling-based method. Graphical methods, in general, mainly focus on the representation of sensitivity analysis results in the form of graphs, charts, or surfaces.

They provide a visual understanding of how an output variable is affected by the variation in input factors. Graphical methods can be used to complement the analysis results of mathematical and statistical methods in establishing systematical framework of sensitivity analysis. Some other classifications of sensitivity analysis methods may focus on the capability of a specific technique, which aids in understanding the applicability of a specific method to a particular model and analysis objective [START_REF] Saltelli | Sensitivity analysis[END_REF].

One-at-a-time method

One-at-a-time method for sensitivity analysis, also known as nominal range sensitivity analysis, local sensitivity analysis or threshold sensitivity analysis, is one of the simplest and most widely used approaches. This method is individually varying only one of the model inputs across its entire range of plausible values at a time while fixing the others at their base-case or mean values, to see what effect this exerts on the model output. The sensitivity of the model output with respect to a particular input variable can be identified by the difference in the model output contributed by the variation of the variable. Regarding deterministic models, sensitivity measures can be typically represented as a positive or negative percentage change of the output variable compared to its normative value.

For linear models, one-at-a-time analysis is advantageous for recognising the most important factors because of its simplicity and low computational cost in implementation.

This approach, however, does not consider the simultaneous variation of input variables.

This makes it fail in identifying the impacts of interaction effects among multiple input factors on the uncertainty in an output variable. Accordingly, the analysis results of one-at-a-time strategy are potentially misleading for models other than linear ones.

One-at-a-time method is quite often used for sensitivity analysis of models with a large number of input variables (e.g., climate models, contagious disease spreading models) since it is easily operated, by only repeating sensitivity analysis process for any number of individual model inputs [START_REF] Happe | Agent-based modelling and sensitivity analysis by experimental design and metamodelling: an application to modelling regional structural change[END_REF][START_REF] Bailis | Mortality and greenhouse gas impacts of biomass and petroleum energy futures in africa[END_REF][START_REF] Campbell | Photosynthetic control of atmospheric carbonyl sulfide during the growing season[END_REF][START_REF] Lopez | Global and regional burden of disease and risk factors, 2001: systematic analysis of population health data[END_REF]. Murphy et al, for example, applied the method to address the range of climate changes resulted from the variation of input factors [START_REF] Murphy | Quantification of modelling uncertainties in a large ensemble of climate change simulations[END_REF].

With one-at-a-time strategy, they obtained statistical predictions of climate prediction and sensitivity indices in assuming that the effects of perturbations of 32 input parameters combine linearly and independently. In a deterministic HIV model with 20 uncertain input parameters involved and supposed to be independent of each other, Blower and Dowlatabadi employed PRCC to assess the statistical relationship between each input parameter and the outcome variable while fixing the remaining input factors at their nominal values [3].

Regarding one-at-a-time method for sensitivity analysis, one often considers the Morris method [START_REF] Morris | Factorial sampling plans for preliminary computational experiments[END_REF]. Morris method is stated as effective to screen a few important input factors from a large number controlling a model. In this method, the input space (for simplicity, defined as an n-dimensional unit hypercube) is discretized to n-dimensional p-level grid.

Each input x i may take values from a sequence {0, 1/(p -1), 2/(p -1), • • • , 1}. For a given value of input vector x, the elementary effect of the ith input factor is defined based on Morris method as

d i (x) = f (x + e i ∆) -f (x) ∆ = ∂f ∂x i , (1.4)
where ∆ is a predetermined multiple of 1/(p -1), and e i a column vector where the ith entry is set to 1 and the rest ones are set to zero. The finite distribution of elementary effects is estimated by randomly sampling different x from input space. For input x i , the distribution is denoted by F i . The mean value u i and standard deviation υ i of the distribution F i are then estimated. u i characterises the effect of input x i on the model output and υ i the nonlinear effect of x i as well as the interaction effects associated with

x i :

u i = ∫ ∂f ∂x i dx, υ i = [ ∫ ( ∂f ∂x i -u i ) 2 dx ] 1/2 . (1.5)
When the model under discussion is non-monotonic, distribution F i contains positive and negative elements. Averaging rule may cancel some effects so as to make u very small or even zero. For this reason, an improved sensitivity measure is considered by Campolongo et al, called u * , which is defined as the mean of the distribution of the absolute values of the elementary effects [START_REF] Campolongo | Sensitivity analysis: the morris method versus the variance based measures[END_REF][START_REF] Campolongo | Enhancing the morris method[END_REF]:

u * = ∫ ∂f ∂x i dx. (1.6)
Measure u * can help identify out the input factors of important overall influence on the model output.

Regression analysis method

Regression analysis is a statistical process for providing an algebraic representation of the relationship between output variable and one or more of input parameters [START_REF] Darlington | Regression analysis and linear models: Concepts, applications, and implementation[END_REF]. It allows one to understand how the output variable changes when any one of the input factors is varied while the remaining factors are fixed.

The earliest form of the regression was the method of least squares, which was considered by Legendre in 1805 [START_REF] Legendre | Nouvelles méthodes pour la détermination des orbites des comètes[END_REF] and also by Gauss in 1809 [START_REF] Gauss | Theory of the motion of the heavenly bodies moving about the sun in conic sections: a translation of Carl Frdr[END_REF]. However, the term "regression" was proposed by Francis Galton in the late of nineteenth century, with which, a biological system was described [START_REF] Galton | Typical laws of heredity[END_REF]. It was later extended by Udny Yule and Karl Pearson to a more general statistical context [START_REF] Yule | On the theory of correlation[END_REF].

In the context of sensitivity analysis, regression analysis usually involves the construction of linear relationship connecting output variable and input parameters. The standardized regression coefficients are then directly used for assessing the sensitivity of model output with respect to individual input factors. Regression analysis contains three groups of variables: the unknown regression coefficients, denoted as

b i with i = 0, 1, • • • , n, the input factors x i with i = 1, 2, • • • , n
, and the model output y. y could be a vector.

But for simplicity, we consider y as a scalar.

Regression analysis is most properly performed by independent random samples which constitute the mapping from input factors to the output variable. By linear regression, the model under study is approximated as

ŷ = b 0 + n ∑ i=1 b i x i , (1.7)
where b i is the regression coefficient for input x i , which can be interpreted as the change in output y when the input factor x i increases by one unit in keeping the remaining factors constant [START_REF] Peck | Statistics: The exploration & analysis of data[END_REF]; ŷ denotes the predicted value of output variable for a given point in the n-dimensional input space when regression coefficients are determined. The coefficients b i are determined by least squares: minimizing the sum of squares of deviation from the true values:

M ∑ j=1 (y j -ŷj ) 2 = M ∑ j=1 [ y j - ( b 0 + n ∑ i=1 b i x i j )] 2 , (1.8)
where M is the number of samples (experimental points), y j the jth output data point given by the jth n-dimensional input data point, x i j the jth sampled value of input x i [START_REF] Freedman | Statistical models: theory and practice[END_REF]. The deviation of the prediction of the regression model (Eq. (1.7)) from the exact values given by the original model can be evaluated as

R 2 = M ∑ j=1 (ŷ j -ȳ) 2 / M ∑ j=1 (y j -ȳ) 2 , (1.9)
which provides a measure of the amount of uncertainty in output variable explained by linear regression model [START_REF] Draper | Applied regression analysis[END_REF][START_REF] Montgomery | Introduction to linear regression analysis[END_REF]. Particularly, R 2 → 1 indicates that the developed regression model accounts for most of the uncertainty in output variable. Conversely, R 2 → 0 means that the regression model is not satisfied in explaining the output uncertainty [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF].

To some degree, the regression coefficients can reflect the sensitivity of model output to input factors. If a coefficient b i is close to 0, then there is not a statistically significant linear relationship between input x i and the output y. Conversely, if b i is significantly different from 0, then the output y can be regarded as being sensitive to x i . However, b i is influenced by the units of x i . To reduce the dimensional effects of input factors, the regression model represented by Eq. (1.7) is commonly standardised to make the variance of output and input variables equal to 1:

(ŷ -ȳ)/σ = n ∑ i=1 (b i σi /σ)(x i -µ i )σ i , (1.10) where σ
=   1 M M ∑ j=1 (y j -ȳ) 2   1/2 , (1.11) σi =   1 M M ∑ j=1 (x i j -µ i ) 2   1/2 , (1.12) and ȳ
= 1 M M ∑ j=1 y j , (1.13)
µ i = 1 M M ∑ j=1 x i j . (1.14)
The coefficients b i σi /σ are referred to as standardised coefficients, taking values between -1 and 1. The standardised coefficients are helpful in identifying which of input parameters have greater effects on the output variable when the input variables are measured in different units of measurement [START_REF] Neter | Applied linear statistical models, volume 1 of Irwin series in statistics[END_REF]. Linear regression analysis is most suitable when the model of interest is in fact linear as it is difficult to interpret the standardised coefficients when nonlinear regression analysis is involved. In some analyses, nonlinear regression provides an alternative to linear regression for more accurate estimation of the relationship between output and input variables [START_REF] Bates | Nonlinear regression analysis and its applications[END_REF][START_REF] Motulsky | Fitting models to biological data using linear and nonlinear regression: a practical guide to curve fitting[END_REF].

Because of simplicity and low computational cost, regression analysis as a strategy of sensitivity analysis has been adopted by many researchers from various of fields, such as medical science [START_REF] Prentice | Regression analysis of grouped survival data with application to breast cancer data[END_REF][START_REF] Taylor | Alteration in lymphocyte subpopulations in bovine leukosis virus-infected cattle[END_REF], bioscience [START_REF] Daniel | Biostatistics: a foundation for analysis in the health sciences[END_REF], human science [START_REF] Happe | Agent-based analysis of agricultural policies: an illustration of the agricultural policy simulator agripolis, its adaptation and behavior[END_REF], and food science [START_REF] Torjusen | Food system orientation and quality perception among consumers and producers of organic food in hedmark county, norway[END_REF][START_REF] Helton | Uncertainty and sensitivity analysis of food pathway results with the maccs reactor accident consequence model[END_REF].

Response surface method

Response surface method consists of a group of mathematical and statistical techniques used in the development of an adequate model function connecting an output variable and a number of input parameters. With the established functional relationship, response surface method can identify curvature in the response surface by accounting for high-order effects produced by input parameters. The method was introduced by Box and Wilson in 1951 [START_REF] Box | On the experimental attainment of optimum conditions[END_REF]. The main idea of this approach is to use a sequence of designed experiments to obtain an optimal response. Considering the complexity in implementation, response surface method, therefore, is commonly used for the analysis of models with limited number of input factors.

In general, the functional relationship between model output and input parameters is unknown but can be approximated by a low-degree polynomial model of the form

y = f T (x)β + ϵ, (1.15)
where

x = (x 1 , x 2 , • • • , x n ) T , the input vector of n-dimensional variables; f T (x) is a
vector function of a group of elements, consisting of powers and cross-products of powers of individual input parameters up to a certain degree d (≥ 1); β is a vector of unknown constant coefficients; ϵ is a random experimental error and assumed to have a zero mean [START_REF] Khuri | Response surface methodology[END_REF].

Currently two important models are used in response surface method, with one being linear and the other, nonlinear. The linear one is also classified as the first-degree model (d=1) constructed in terms of the first-order terms of input parameters:

y = β 0 + n ∑ i=1 β i x i + ϵ. (1.16)
The nonlinear one is classified as the second-degree model (d=2), which still involves second-order effects of input parameters, except for first-order ones:

y = β 0 + n ∑ i=1 β i x i + n ∑ i n ∑ j>i β ij x i x j + n ∑ i=1 β ii x 2 i + ϵ.
(1.17)

The application of response surface method to sensitivity analysis of models can be concluded as three procedures:

(1) To approximately establish the functional relationship connecting output variable and input factors.

(2) To quantify, through hypothesis testing, importance of individual factors.

( In the estimation of the first-degree model, an easy but most common design is 2 n factorial design [START_REF] Box | Statistics for experimenters: design, innovation, and discovery[END_REF]. In a 2 n factorial design, each input variable is measured at two levels which are commonly coded as -1 for the low level and +1 for the high level. A factorial design consists of all possible combinations of previously defined levels of n input factors. In practical applications, the points in a two-level factorial design are frequently represented by plus and minus signs, conventionally, -for the first (or low) level, and + for the second (or high) level. Take the case of three input factors as an example. The corresponding 2 3 design is a 8 × 3 matrix of the form

D =     ---+ -+ + + --+ -+ -+ + -+ --+ + -+     T . (1.18)
If n is large, a large number of points will be introduced by the 2 n factorial design, thereby high cost is required in the computer simulation. For this case, fractions (e.g., one-half fraction, one-fourth fraction) of a 2 n design are often considered to reduce the cost of computer simulations in the estimation of the first-degree model. In general, a 2 -m th fraction of a 2 n design contains 2 n-m points of a 2 n design. Here m is an integer number such that 2 n-m ≥ n+1 for guaranteeing all n+1 parameters (elements of vector β) included in the first-degree model (Eq. (1.16)) can be estimated. Particular manners for the construction of fractions of a 2 n design can be found in Refs. [START_REF] Box | Statistics for experimenters: an introduction to design, data analysis, and model building[END_REF][START_REF] Montgomery | Design and analysis of experiments[END_REF]. Some other commonly discussed designs for fitting the first-degree model are Plackett-Burman design [START_REF] Plackett | The design of optimum multifactorial experiments[END_REF][START_REF] Ahuja | Application of plackett-burman design and response surface methodology to achieve exponential growth for aggregated shipworm bacterium[END_REF] and simplex design [START_REF] Khuri | Response surfaces: designs and analyses[END_REF]. The Plackett-Burman design allows two levels of each input factor, analogous to the 2 n factorial design, but requires a much smaller number of design points, especially for large n. The number of design points required by the Plackett-Burman design is equal to the number of parameters (=n + 1) in the firstdegree model. Specifically, this design can be employed only when the number of input variables, n, is a multiple of 4 [START_REF] Plackett | The design of optimum multifactorial experiments[END_REF]. The simplex design also contains n + 1 experimental points. These points are located at the nodes of an n-dimensional regular-sided figure [START_REF] Box | Multi-factor designs of first order[END_REF][START_REF] Rajashekhar | A new look at the response surface approach for reliability analysis[END_REF].

In the fitting of the second-degree model, one of the most frequently discussed designs is the 3 n factorial design [START_REF] Box | Some new three level designs for the study of quantitative variables[END_REF]. A 3 n factorial design is formed from all possible combinations of the levels of all n input variables. Each input variable has three levels that are commonly coded as -1 for the low level, 0 for the intermediate level, and +1 for the high level. In practice, the matrix of a 3 n factorial design simply consists of plus and minus signs, and also 0. Analogously, the number of experimental points for this design (=3 n ) will be very large when a large number of input factors are involved into the original model. Following the phenomenon, fractions of a 3 n factorial design are often employed to save the cost of computer simulations. In the construction of fractions, the number of experimental points must at leat equal the number of parameters (=2n+1+n(n-1)/2, the number of elements of vector β) included in the second-degree model (Eq. (1.17)) [START_REF] Montgomery | Design and analysis of experiments[END_REF][START_REF] Mclean | Applied factorial and fractional designs. Statistics: a Series of Textbooks and Monogrphs[END_REF]. Another most widely used second-order design is the central composite design. It is stated as the most popular design for building a second-degree model. The central composite design was first introduced by A. I. Khuri in 1988 [START_REF] Khuri | A measure of rotatability for response-surface designs[END_REF], consisting of three distinct sets of experimental runs:

(1) A full (or a fraction of) 2 n factorial design. This is called the factorial portion. Two levels of each input factor are coded as -1 and +1. They are often simplified as plus and minus signs in design matrix.

(2) n 0 central points. Central point, commonly coded as 0, is the median of the values of each factor used in the factorial portion. n 0 replications of central point is used to improve the precision of the experiment.

(3) 2n axial points. Two points are taken on the axis of each input variable at a distance of α from the center of the variable.

The number of experimental runs (design points) in a central composite design is (2 n + 2n + n 0 ). The design matrix of a simple case with n = 3 and n 0 = 4 is formed as follows:

D =     ---+ -+ + + -α α 0 0 0 0 0 0 0 0 --+ -+ -+ + 0 0 -α α 0 0 0 0 0 0 -+ --+ + -+ 0 0 0 0 -α α 0 0 0 0     T . (1.19)
Many strategies have been developed to select a useful value of the axial parameter α.

Let F denote the number of points in the factorial portion and T = 2n + n 0 . Two common values are

α = (Q × F/4) 1/4 (1.20) with Q = ( √ F + T - √ F ) 2
, which makes the central composite design orthogonal, and

α = F 1/4 (1.21)
which makes the design rotatable [START_REF] Myers | Response surface methodology[END_REF]. The value of n 0 is often assigned in terms of some certain desirable properties of the central composite design. For example, n 0 can be set to a value that makes a rotatable central composite design hold orthogonality property or the uniform precision property [START_REF] Khuri | Response surfaces: designs and analyses[END_REF][START_REF] Box | Multi-factor experimental designs for exploring response surfaces[END_REF]. The other most frequently used secondorder design is the Box-Behnken design devised by G. E. P. Box and D. Behnken in 1960 [START_REF] Box | Some new three level designs for the study of quantitative variables[END_REF]. In this design, three levels (equally spaced) of each input factor are considered.

A Box-Behnken design is formed from a particular subset of the full 3 n factorial design [START_REF] Box | Response surfaces, mixtures, and ridge analyses[END_REF][START_REF] Myers | Response surface methodology: process and product optimization using designed experiments[END_REF]. Some other second-order designs are available in Refs. [START_REF] Doehlert | Uniform shell designs[END_REF][START_REF] Box | Factorial designs, the -x'x-criterion, and some related matters[END_REF][START_REF] Hoke | Economical second-order designs based on irregular fractions of the 3 n factorial[END_REF][START_REF] Roquemore | Hybrid designs for quadratic response surfaces[END_REF].

A representative strategy to develop the expression of Eq. (1.15) is using a least squares repression method to fit a standardized first-or second-order equation to the data obtained from the original model. MC methods are typically borrowed to produce multiple values of each input factor, thereby to calculate corresponding values of model output.

Other techniques such as rank-based or nonparametric approaches are also occasionally considered in employing the response surface method to establish the framework of sensitivity analysis [START_REF] Khuri | Response surfaces: designs and analyses[END_REF][START_REF] Vidmar | A monte carlo study of robust and least squares response surface methods[END_REF]. The precision and accuracy of analysis results provided by the response surface method can be evaluated by comparing the predictions provided by the method to the output values of the original model generated by the same values of input parameters. If the precision and accuracy are not satisfactory, an improved fit might be obtained by iterating on values of parameters [START_REF] Gardiner | Experimental design techniques in statistical practice: A practical software-based approach[END_REF].

Time-consuming and effort requirement in applying response surface approach typically rely on the number of input parameters included and the type of response structure required. Hence, mainly focusing on the effects of those input factors that have been identified as quite important through a screening sensitivity analysis method (e.g., oneat-a-time method) may be advantageous to reduce the complexity and difficulty in the implementation of response surface strategy.

A key advantage of the response surface method is that one can save computational time in computationally intensive model run by simplifying the form of the model under discussion. Furthermore, the functional form of the model that is established by the response surface method and the values of the coefficients included in the form can provide a fruitful information for quantifying the sensitivity of model output with respect to individual parameters. However, most frameworks established by the response surface method only consider the effects of some but not all of the input factors contained in the original model. This may result in absent or non-accurate global sensitivity measures in sensitivity analysis.

Differential-based method

Differential techniques for sensitivity analysis, also referred to as the direct or local methods, involve partial derivatives of output variable with respect to input parameters. In sensitivity analysis, one of the most used differential-based strategies is the first-order Taylor series approximation. It was discussed since 1966 by engineers and experimentalists [START_REF] Ku | Notes on the use of propagation of error formulas[END_REF].

Recall a generic model of the form y = f (x) with x = (x 1 , x 2 , • • • , x n ) T labeling the input vector of n-dimensional variables. By employing the first-order Taylor series approximation, the variance of output y, denoted as V (y), is calculated as

V (y) = n ∑ i=1 ( ∂y ∂x i ) 2 X 0 V (x i ), (1.22)
where the subscript X 0 indicates that the derivative is taken at a fixed point (often indicated by the central point) in the space of input variables, and V (x i ) the variance of input x i . The sensitivity coefficient, denoted by s i , interprets the importance of input

x i in establishing the uncertainty of output y. s i is determined by (see [START_REF] Helton | Sensitivity analysis of the asymptotic behavior of a model for the environmental movement of radionuclides[END_REF])

s i = ( ∂y ∂x i ) 2
This calculation is performed under the assumption that high-order (≥ 2) partial differentials are negligible and input parameters are independent of each other. Consequently, the first-order Taylor series approximation can provide accurate and reliable analysis results only when the model under study is almost linear or the uncertainty in input parameters is negligible [START_REF] Padulo | Comparative analysis of uncertainty propagation methods for robust engineering design[END_REF][START_REF] Koda | Sensitivity analysis of partial differential equations with application to reaction and diffusion processes[END_REF].

Sensitivity analysis proposed by the differential-based technique is computationally efficient but bound with intensive effort requirement in solving differential equations. When an explicit algebraic equation describes the relationship connecting model output and input parameters, it is straightforward to evaluate sensitivity measures by the use of differential-based strategy. If a large set of equations are involved with the model under discussion, the first-order partial derivative can be approximated as a finite variation in output values driven from a small change in the input parameter [START_REF] Downing | An examination of responsesurface methodologies for uncertainty analysis in assessment models[END_REF][START_REF] Hamby | A review of techniques for parameter sensitivity analysis of environmental models[END_REF]. By neglecting non-linearities of models, the sensitivity of model output with respect to an arbitrary input parameter x i can be approximated as

s i = %∆y %∆x i = [f (x + ∆ i ) -f (x)] /f (x) [x i -∆ i ] /x i , (1.24)
where

(x+∆ i ) = (x 1 , • • • , x i +∆ i , • • • , x n )
, and ∆ i is a small change introduced to input

x i .
A derivative-based global sensitivity method has also been proposed by Sobol and Kucherenko, by averaging the square of local derivatives [START_REF] Sobol | Derivative based global sensitivity measures and their link with global sensitivity indices[END_REF][START_REF] Sobol | Derivative based global sensitivity measures[END_REF]. In this method, the global sensitivity measures are defined as

v i = ∫ C n ( ∂f ∂x i ) 2 dx, (1.25)
where

C n = (x|0 ≤ x i ≤ 1; i = 1, 2, • • • , n), the n-dimensional unite hypercube. v i can
be regarded as an improvement of the importance criterion u * (see Eq. (1.6)). The above definition is motivated by the fact that a high value of the derivative of model output with respect to an input variable indicates a robust influence of the input variable on the model output [START_REF] Lamboni | Derivative-based global sensitivity measures: general links with sobolindices and numerical tests[END_REF]. It is proved that

s T i ≤ v i π 2 V (y) , (1.26)
where s T i are the one-dimensional total sensitivity indices (see Eq. (1.38)) and V (y) the total variance of output y. This states that small v i imply small s T i . Unessential input parameters then can be identified out based on computed values of

v i (i = 1, 2, • • • , n).
For highly nonlinear functions, however, the ranking of influential parameters in terms of the importance criterion v i may suggest false conclusions [START_REF] Sobol | Derivative based global sensitivity measures and their link with global sensitivity indices[END_REF].

The differential-based strategy is usually more demanding than other methods in the sensitivity analysis of complex models. It requires of model designers to explicitly calculate the first-order partial derivatives of output variable with respect to individual parameters, and yet provides only comparable but not accurate results, especially for nonlinear models.

Variance-based methods

Variance-based techniques have a long history in the aspect of sensitivity analysis. They are often used for determining whether an output variable is statistically associated with one or more input factors, and whether the values of model output vary in a statistically significant manner with the variation in values of one or more input variables. In the seventies, Cukier firstly established variance-based sensitivity analysis of multi-variate systems by Fourier implementation [START_REF] Cukier | Nonlinear sensitivity analysis of multiparameter model systems[END_REF]. While the complete variance decomposition strategy was firstly developed by Sobol in 1993 [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. In 1994, Jansen et al introduced an efficient method relying on random sampling to evaluate the partial contributions from input variables of uncertainty to the predicted variance in output variable [START_REF] Jansen | Monte carlo estimation of uncertainty contributions from several independent multivariate sources[END_REF]. A similar strategy was developed by Homma and Saltelli [START_REF] Homma | Importance measures in global sensitivity analysis of nonlinear models[END_REF] in 1996 to determine global sensitivity measures that quantify the global importance of individual input variables in the estimation of model response. In spite of time consuming in computation, the instrument of complete variance decomposition is known to be useful and informative for uncertainty and sensitivity analysis of complex nonlinear systems [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF].

Variance-based sensitivity analysis, often referred to as the Sobol method or Sobol indices, specifies the uncertainty in input and output variables through probability distributions. Working within a probabilistic framework, it decomposes the variance of model output into different partial contributions attributable to individual input variables and to their combinations. By computing the percentage of each partial variance contribution in the global variance of output variable, sensitivity measures are directly interpreted for individual factors and also for their interaction effects. Variance-based sensitivity analysis methods are attractive and widely used because they allow full exploration of input space, analysis of nonlinear models, and consideration of interactions between different input variables.

Without any assumption regarding the type of the model under discussion, variancebased approaches find broad applications across various fields, including scientific models evaluation [START_REF] Saltelli | Sensitivity analysis in practice: a guide to assessing scientific models[END_REF], risk assessment [START_REF] Gerken | Models for retrospective quantification of indoor radon exposure in case-control studies[END_REF], importance assessment [START_REF] Saltelli | Sensitivity analysis for importance assessment[END_REF], economic system analysis [START_REF] Zhu | An analytic method for sensitivity analysis of complex systems[END_REF], behaviour prediction in forest systems [START_REF] Song | Variance-based sensitivity analysis of a forest growth model[END_REF], etc.

Recall the generic model of the form y = f (x) with x = (x 1 , x 2 , • • • , x n ) T labeling the input vector of n-dimensional variables. The model is defined over C n , the n-dimensional unit hypercube, as defined before. Recalling the classical Hoeffding decomposition [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF][START_REF] Efron | The jackknife estimate of variance[END_REF][START_REF] Sacks | Design and analysis of computer experiments[END_REF][START_REF] Archer | Sensitivity measures, anova-like techniques and the use of bootstrap[END_REF][START_REF] Rabitz | Efficient inputoutput model representations[END_REF], the output variable can be expanded as

y = f 0 + n ∑ i=1 f i (x i ) + n ∑ i=1 n ∑ j>i f ij (x i , x j ) + • • • + f 12•••n (x 1 , x 2 , • • • , x n ), (1.27)
where f 0 is a constant, f i a function of x i , f ij a function of x i and x j , and so on up to the last term a function involving all input variables. Each term is square integrable over C n . Summands presented in Eq. (1.27) must satisfy

∫ 1 0 f i 1 i 2 •••is (x i 1 , x i 2 , • • • , x is )dx i k = 0, (1.28)
where

1 ≤ i 1 < i 2 < • • • < i s ≤ n and i k ∈ {i 1 , i 2 , • • • , i s }. This condition drives f 0 = ∫ C n f (x)dx, (1.29)
f i (x i ) = ∫ 1 0 • • • ∫ 1 0 f (x)dx/dx i -f 0 , (1.30)
f ij (x i , x j ) = ∫ 1 0 • • • ∫ 1 0 f (x)dx/dx i dx j -f 0 -f i (x i ) -f j (x j ), (1.31)
. . . where dx/dx i is the product of all the dx j except dx i , similar to dx/dx i dx j which indicates the integration with respect to all input variables, except x i and x j . By assuming f (x) is square-integrable, the variance of model output can be expressed as

V (y) = ∫ 1 0 • • • ∫ 1 0 [ n ∑ s=1 n ∑ i 1 <•••<is f 2 i 1 ,••• ,is (x i 1 , x i 2 , • • • , x is )dx i 1 • • • x is ] . (1.32)
Expanding the above equation yields

V (y) = n ∑ i=1 V i + n ∑ i=1 n ∑ j>i V ij + • • • + V 12•••n , (1.33)
where

V i = ∫ 1 0 f 2 i (x i )dx i , (1.34)
V ij = ∫ 1 0 ∫ 1 0 f 2 ij (x i , x j )dx i dx j , (1.35)
. . . are in sequence explaining the variances contributed by x i alone, by the interaction between x i and x j , etc. Equation (1.33) is known as the complete variance decomposition, deduced by Sobol in 1993. It shows how the variance of a model output can be decomposed into terms that are attributable to individual input factors and to their interaction effects [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. Following the concept of Sobol's variance decomposition, sensitivity measures are defined by

s i = V i / V (y), (1.36)
s ij = V ij / V (y), (1.37)
. . .

s T i =   V i + n ∑ j=1;j̸ =i V ij + • • • + V 12•••n   / V (y), (1.38)
where s i is often called the first-order sensitivity index or the main effect index, labeling the fraction of V (y) contributed by x i alone; s ij the second-order sensitivity index, indicating the fraction of V (y) contributed by the interaction between x i and x j ; and so on up to the last label s T i the total sensitivity index, quantifying the fraction of V (y) contributed by x i alone and also by interactions of x i with the remaining factors. They

satisfy n ∑ i=1 s i + n ∑ i=1 n ∑ j>i s ij + • • • + s 12•••n = 1, (1.39) n ∑ i=1 s T i ≥ 1. (1.40)
The equal sign in Eq. (1.40) holds iff the model under analysis is purely additive.

Theoretically, the determination of sensitivity measures depends upon multidimensional integrals. In practice, sampling-based strategies are often employed to carry out this mission. It can become computationally expensive when there are a large number of input variables.

Sobol's definitions of the second and higher order partial variance contributions in Eq.

(1.33) hold if input factors are independent of each other. In the presence of input correlations, however, partial variance contributions with dimensionality larger than 1 are contributed not only by the coupling items presented in the functional form of the model under discussion (for independent case), but also by the input correlations.

Recently, the interest in extending sensitivity analysis strategies from uncorrelated case to the correlated one is increasing as correlated input variables are of frequent occurrence in practical applications [START_REF] Jacques | Sensitivity analysis in presence of model uncertainty and correlated inputs[END_REF][START_REF] Veiga | Local polynomial estimation for sensitivity analysis on models with correlated inputs[END_REF]. In general, variance-based sensitivity analysis in the presence of correlated inputs is carried out by using linear correlation model which regards the correlation part of an arbitrary variable as a linear combination of the rest variables [START_REF] Most | Variance-based sensitivity analysis in the presence of correlated input variables[END_REF][START_REF] Hao | A new interpretation and validation of variance based importance measures for models with correlated inputs[END_REF][START_REF] Li | General validation and decomposition of the variance-based measures for models with correlated inputs[END_REF][START_REF] Zhang | A new framework of variance based global sensitivity analysis for models with correlated inputs[END_REF].

Moment independent method

The moment independent method looks at the influence of input uncertainty on the entire output distribution without reference to a specific moment of the output. The moment independent importance indicator was first introduced by Chun et al [START_REF] Chun | An uncertainty importance measure using a distance metric for the change in a cumulative distribution function[END_REF].

They quantified the entire change of CDFs in terms of the normalised Euclidean metric distance between two CDFs. The metric distance (Minkowski distance) of order a be-

tween two points X 1 = (x 1 1 , x 2 1 , • • • , x n 1 ) and X 2 = (x 1 2 , x 2 2 , • • • , x n 2 ) is in general defined by D = ( n ∑ i=1 |x i 1 -x i 2 | a ) 1/a , (1.41)
where a is a number no less than 1. The metric distance is typically used with a being 1 or 2, which correspond to the Manhattan distance and the Euclidean distance, respectively.

The measure of uncertainty importance is then defined by

MD(i : o) = (∫ 1 0 [Y i j -Y o j ] 2 dj ) 1/2 / E(y o ), (1.42)
where Y o j is the jth percentile of the output CDF for the base case, Y i j the jthe percentile of the output CDF after introducing a change to input x i , and E(y o ) the mean value of output y for the base case. The base case refers to the case where an output distribution is obtained by setting all input distributions to their nominal ones. It is stated that a larger value of MD(i : o) implies a more important parameter x i , compared to other input parameters. However, the value of MD(i : o) depends on the hypothesised change introduced to parameter x i .

Another moment independent importance measure was presented by Borgonovo, which does not require one pre-suppose any changes associated with input parameters [START_REF] Borgonovo | Global sensitivity analysis in inventory management[END_REF].

It follows from his concept that the global sensitivity index of parameter x i , denoted by δ i , with respect to the output y is represented by

δ i = 1 2 ∫ P (x i ) [∫ |P y (y) -P y|x i (y)|dy ] dx i , (1.43)
where

P (x i ) = ∫ • • • ∫ P (x)
∏ j̸ =i dx j the marginal density (PDF) of input x i , P (y) is the PDF of output y, and P Y |x i (y) indicates the conditional PDF of y by assuming that the input parameter x i is fixed at a constant value. It is proved that δ i holds some properties, such as a) 0 ≤ δ i ≤ 1; b) δ i = 0 if the output y is independent of x i ; c) δ 1,2,••• ,n = 1 (the importance of all input parameters equals one); etc [START_REF] Borgonovo | A new uncertainty importance measure[END_REF]. The indicator δ helps identify out input parameters that affect output uncertainty the most. In identifying the less relevant parameters with respect to the model output, δ agrees with Sobol's variancebased method. However, discrepancy between two evaluations exists in ranking relevant parameters.

Regarding the difficulty in deriving the PDFS of output variable, a new moment independent measure was recently proposed, called the PAWN index, with the conditional and unconditional distributions characterised by their CDFs [START_REF] Pianosi | A simple and efficient method for global sensitivity analysis based on cumulative distribution functions[END_REF]. The Kolmogorove-Smirnov statistic ( [START_REF] Kolmogorov | Sulla determinazione empirica di una leggi di distribuzione[END_REF]) is employed to measure the distance between conditional and unconditional CDFs:

KS(x i ) = max y |F y (y) -F y|x i (y)|, (1.44)
where F y (y) is the unconditional CDF of output y, F y|x i (y) the conditional CDF when the input x i is fixed. As KS(x i ) depends on the value at which x i is fixed, the PAWN index, denoted by T i , is then defined by considering a statistic (e.g. maximum or median) over all possible values of x i :

T i = stat x i {KS(x i )}. (1.45)
The PAWN index T i is a global, quantitative and model independent sensitivity index, varying between 0 and 1 (the higher the value, the more influential x i ). T i

(i = 1, 2, • • • , n
) can be used for ranking input parameters according to their contributions to the output uncertainty, so as to identify which of input parameters are influential, and which are non-influential. Compared to other moment independent sensitivity indices, the PAWN index has the advantage of being very easy to implement and interpret [START_REF] Pianosi | A simple and efficient method for global sensitivity analysis based on cumulative distribution functions[END_REF]. Regarding the convergence rate and screening of non-influential parameters, the PAWN method is in good agreement with Sobol's variance-based method.

However, the PAWN method is better than Sobol's one in quantifying the difference between relative importance of influential parameters [START_REF] Zadeh | Comparison of variance-based and moment-independent global sensitivity analysis approaches by application to the swat model[END_REF].

Sampling-based method

In the implementation of uncertainty and sensitivity analysis, the sampling-based (also known as MC) approach has been classified as both effective and widely used. In the absence of input correlations. The second one is the generation of sample points in n-dimensional input space according to the previously defined distributions of individual input variables. Diversified sampling procedures have been designed. In this section, we discuss some frequently considered sampling strategies: random sampling, sampling for iterated fractional factorial design (IFFD), importance sampling, Latin Hypercube sampling (LHS), and Quasi Monte Carlo (QMC) sampling.

Random sampling

Sample points provided by random sampling are generated by pseudo-random numbers.

Values of single uncertain input variables are chosen randomly and entirely by chance from their real ranges, independent of their PDFs [START_REF] Starnes | The practice of statistics[END_REF]. It is stated that random sampling is prone to clustering that yields wanted sample points and empty areas in input space because each new sample is selected randomly without minding the gaps between already generated sample values [START_REF] Mckay | Comparison of three methods for selecting values of input variables in the analysis of output from a computer code[END_REF][START_REF] Macdonald | Comparison of sampling techniques on the performance of montecarlo based sensitivity analysis[END_REF].

IFFD

The factorial sampling and fractional factorial sampling have been briefly mentioned before. They are widely considered in the response surface method. Sampling for IFFD was also proposed as a direct tool for sensitivity analysis [START_REF] Andres | Statistical sampling strategies[END_REF][START_REF] Andres | Using iterated fractional factorial design to screen parameters in sensitivity analysis of a probabilistic risk assessment model[END_REF]. In this sampling, input variables are sampled at discrete levels. This is the main difference from other sampling strategies where each input variable is sampled at a continuous interval. IFFD was built to identify a few robust parameters within batches of a large number of fragile ones [START_REF] Andres | Using iterated fractional factorial design to screen parameters in sensitivity analysis of a probabilistic risk assessment model[END_REF]. Investigation indicates that the sensitivity measure based on IFFD is extremely reproducible, and more robust than rank regression coefficient in that it can also identify quadratic effects of individual parameters, except for linear ones [START_REF] Saltelli | Sensitivity analysis of model output. performance of the iterated fractional factorial design method[END_REF].

Specifically, IFFD is generated by multiple iterations of a basic orthogonal fractional factorial design. Given a basic fractional factorial design, three steps of randomization control the generation of a sampling for IFFD:

(1) Parameters are randomly assigned to columns of the basic fractional factorial design.

(2) Each parameter is randomly oriented in considering two orientations: positive one and negative one. Parameters with a positive orientation will copy their values from the associated column of the basic design, and those with a negative one will take the opposite values from the associated column of the basic design.

(3) The orientation value of each parameter is randomly set to zero in a previously defined proportion of the total iterations.

The first two randomization steps are carried out independently for each iteration, and the third one is performed within the whole set of iterations. Each iteration of an IFFD is analysed separately. The obtained results are then combined together for the sensitivity analysis of the system under discussion. Denote the value of the output variable in the ith simulation of the mth iteration by y m [i], and the input parameter that copies values from the jth column of the basic design by z j . The main effect in the mth iteration of parameter z j on the output variable is given by

M E m (z j , y m ) = 1 n 2n ∑ i=1 J n [i, j] • y m [i], (1.46)
where n is the number of input parameters, J n the basic fractional factorial design of size n × n, and J n [i, j] the (i, j)th position in matrix J n . M E m (z j , y m ) is a linear effect as it calculates the difference in average response between two levels (low and high) of

z j .
The main effect of a parameter x throughout the entire design of M iterations is given by

M E(x, y) = ∑ M m=1 S m x • M E(z c m x , y m ) ∑ M m=1 |S m x | , (1.47)
where S m x is the orientation value of parameter x in the mth iteration, c m x (taking a value from 1 to n) the randomly chosen column associated with parameter x. The quadratic effect of parameter x can also be defined by

QE(x, y) = ∑ M m=1 (1 -|S m x |) ∑ 2n i=1 y m [i] 2n ∑ M m=1 (1 -|S m x |) - ∑ M m=1 |S m x | ∑ 2n i=1 y m [i] 2n ∑ M m=1 |S m x |
.

(1.48)

Importance sampling

In statistics, importance sampling is a well-known technique for estimating properties of a particular distribution, by using samples generated from a different distribution than the distribution of interest. It is often used for the estimate of rare event probability [START_REF] Morio | An overview of importance splitting for rare event simulation[END_REF][START_REF] Morio | Importance sampling: how to approach the optimal density?[END_REF]. Consider a general univariate function f (x). x is a random variable from a previously defined distribution of density P (x). One wishes to estimate the expected value of f (x) under P (x), denoted as E(f ; P ). A MC estimate of E(f ; P ) is

Ê(f ; P ) = ∫ f (x)P (x)dx = 1 M M ∑ j=1 f (x j ) (1.49)
where M denotes the sample size. The precision of this estimate depends upon the variance of f (x):

V ( Ê(f ); P ) = V (f, P )/M. (1.50)
The basic idea of importance sampling is to reduce V ( Ê(f ); P ) for a given sample size by using samples from a different distribution. The principle of importance sampling is

Ê(f ; P ) = ∫ f (x)P (x)dx = ∫ f (x) P (x) G(x) G(x)dx, (1.51)
where the newly introduced distribution G(x) should have the same support as P (x).

Apparently, sampling x from a distribution P (x) is equivalent to sampling x • w(x) from the new distribution G(x), with importance sampling weight:

w(x) ≡ P (x) G(x) . (1.52)
The object in importance sampling is to concentrate the distribution of the sample points in the parts of interval of most importance (e.g. the interval [0, 1] that we discussed in the example) instead of spreading them over the whole region [START_REF] Hammersley | General principles of the monte carlo method[END_REF]. A good importance sampling function G(x) should hold the following properties [START_REF] Anderson | Monte carlo methods and importance sampling[END_REF] (1) G(x) > 0 whenever f (x) ̸ = 0;

(2) G(x) should be most closely the shape of f (x) and makes f (x) • P (x)/G(x) bounded;

(3) It should be easy to simulate values from G(x);

(4) It should be easy to compute the density G(x) for any values of x that one might realize. Consider a simple example of f (x) = x with x following a standard normal distribution:

x ∼ N (0, 1). One wants to evaluate the mean of f in the region [0, 1] of x. By MC sampling, a sequence values

x j (j = 1, • • • , M ) of x can be generated from N (0, 1). The mean value of f (x) is then calculated by Ê(f ; P ) = ∑ M j=1 c j • x j ∑ M j=1 c j , (1.53)
with

c j = { 1 0 ≤ x j ≤ 1, 0 otherwise. (1.54)
In importance sampling, we set

G(x) = 1 for x ∈ [0, 1]
. M sampled values of x then can be obtained from a uniform distribution in the region [0, 1]. For each sampled value x j , the importance weight is .55) where P (•) is the standard normal distribution. The mean value of f (x) based on importance sampling is then calculated by

w(x j ) = P (x j ) ∫ 1 0 P (z)dz . ( 1
Ê(f ; P ) = 1 M M ∑ j=1 x j • w(x j ). (1.56)
A comparison between MC and importance sampling is presented in table 1.1. Results display that the variance of the mean value of f (x) is reduced by introducing importance sampling.

LHS

LHS is a kind of stratified sampling strategy. It is widely regarded as one of the most popular variance reduction techniques that can be employed to increase the analysis efficiency so as to more possibly achieve the desired accuracy. The LHS was firstly adjacent intervals with equal probability provided by P 1 . One point is selected at random from each interval to construct a sequence of M points {x

1 j }, j = 1, 2, • • • , M . Another sequence, {x 2 j }, j = 1, 2, • • • , M
, can be formed similarly but independently for x 2 . The two generated sequences are then randomly paired without replacement to constitute an M × 2 matrix. These M pairs are continuously randomly combined without replacement with M points in the sequence {x 3 j }, j = 1, 2, • • • , M , to form an M × 3 matrix, and so on until an M × n matrix is constructed: The above generated LHS is known as standard LHS (or randomly generated LHS). It can be quite structured: the variables may be highly correlated or the sampled values may not have good uniformity of distribution in the space of input variables. Consequently, some optimal criteria are proposed to avoid these problems and to obtain an optimal sampling which achieves the space-filling property of input factors.

x j = [x 1 j , x 2 j , • • • , x n j ], j = 1, 2, • • • , M. (1.57)
An often considered optimization is designed based on the maximin distance criterion [START_REF] Johnson | Minimax and maximin distance designs[END_REF]. For an established sampling design D (an M × n matrix), the iner-site distance between any two n-dimensional sample points (two rows in the M ×n matrix, as displayed in Eq. (1.57)) is defined by

d(x t , x s ) = [ n ∑ i=1 (x i t -x i s ) α ] 1/α , (1.58)
where α takes two values: 1 (corresponding to the rectangular distance) or 2 (related to the Euclidean distance). The first built maximin distance criterion aims at maximizing the minimum inter-site distance:

min 1≤t,s≤M ;t̸ =s d(x t , x s ). (1.59)
Later on, an intuitively appealing extension of the maximin distance criterion was introduced by Morris and Mitchell [START_REF] Morris | Exploratory designs for computational experiments[END_REF]. For a designed sampling D, they defined a distance list:

(d 1 , d 2 , • • • , d m )
where the elements are the distinct values of iner-site distances and sorted with

d 1 < d 2 < • • • < d m .
Obviously, m, labeling the number of distinct distance values, can be as large as

( M 2
)
(the iner-site distances between any two sampling items are unique) and as small as 1 (the iner-site distances between any two sampling items are consistent). Continuously, an index list was accordingly defined:

(J 1 , J 2 , • • • , J m )
in which an element J k counts the number of pairs of sites in the designed sampling separated by d k . The maximin distance design criterion tends to make the minimum distance maximal. This implies a sampling driven by this criterion must follow the property: (1) d 1 is maximized; (2) J 1 is minimized. To construct an optimal sampling that satisfies this property, an evaluation function was then introduced:

ϕ p (D) = [ m ∑ k=1 J k d -α k ] 1/α , (1.60)
where α is a positive integer and D a designed sampling. An optimal sampling is generated by minimising Eq. (1.60) [START_REF] Morris | Exploratory designs for computational experiments[END_REF].

The Shannon information entropy defined by Shannon in 1948 was also applied to the optimisation of LHS, known as the entropy criterion [START_REF] Shannon | A mathematical theory of communication[END_REF]. Shannon demonstrated that the lower the information entropy, the more precise the knowledge is. It has been further illustrated that the entropy criterion is equivalent to minimizing the quantity:

-log |R|, (1.61)
where R is the correlation matrix of the designed sampling with elements given by [START_REF] Owen | Orthogonal arrays for computer experiments, integration and visualization[END_REF]. The entropy criterion can help design an optimal sampling that minimizes correlations among input factors.

r ts = exp [ n ∑ i=1 c i |x i t -x i s | α ] , 1 ≤ t, s ≤ M ; 1 ≤ α ≤ 2, (1.62) in which c i (i = 1, 2, • • • , n) are correlation coefficients
In experimental design, L 2 discrepancy is most frequently used in solving the nonuniformity problem. The centered L 2 -discrepancy (CL 2 ) criterion was developed to optimise both random sampling and LHS [START_REF] Fang | Uniform design: theory and application[END_REF]. According to this criterion, a constructed sampling D (an M × n matrix) is optimal (having space-filling property) if it minimizes the centered L 2 discrepancy given by

CL 2 (X) = ( 13
) 2 - 2 M M ∑ j=1 n ∏ i=1 ( 1 + 1 2 |x i j -0.5| - 1 2 |x i j -0.5| 2 ) + 1 M 2 M ∑ j=1 M ∑ k=1 n ∏ i=1 ( 1 + 1 2 |x i j -0.5| + 1 2 |x i k -0.5| - 1 2 |x i j -x i k | 12
) .

(1.63) Some other early developed algorithms also exist in constructing optimal LHSs. They have been shown to have a good space-filling property, including, for example, the rowwise element exchange algorithm [START_REF] Park | Optimal latin-hypercube designs for computer experiments[END_REF], the columnwise-pairwise algorithm [START_REF] Kenny | Algorithmic construction of optimal symmetric latin hypercube designs[END_REF], and the threshold accepting algorithm [START_REF] Fang | Centered l2-discrepancy of random sampling and latin hypercube design, and construction of uniform designs[END_REF]. Considering the requirement of highly computational cost in the implementation of these algorithms, Jin et al developed an efficient and flexible method for constructing optimal samplings, by the use of enchanced stochastic evolutionary algorithm and the employment of efficient evaluating strategies of the optimality criteria [START_REF] Jin | An efficient algorithm for constructing optimal design of computer experiments[END_REF]. Some recently introduced methods for optimising the strategy of LHS can be found in Refs. [START_REF] Joseph | Orthogonal-maximin latin hypercube designs[END_REF][START_REF] Viana | An algorithm for fast optimal latin hypercube design of experiments[END_REF][START_REF] Pholdee | An efficient optimum latin hypercube sampling technique based on sequencing optimisation using simulated annealing[END_REF].

QMC sampling

QMC method is specifically designed by using low-discrepancy sequences (also called quasi-random sequences or sub-random sequences) to place sampled values as uniformly as possible. This is in contrast to the regular MC method (random sampling) consisting of sequences of pseudo-random numbers. MC and QMC methods are stated in a similar way. The problem is to approximate the integral of a function as the average of the function evaluated at a set of sampled values. The difference between QMC and MC methods is the way sampled values are generated. There are a few well-known and commonly used low-discrepancy sequences, including Halton sequences [START_REF] Halton | Algorithm 247: Radical-inverse quasi-random point sequence[END_REF], Sobol sequences [START_REF] Joe | Remark on algorithm 659: Implementing sobol's quasirandom sequence generator[END_REF], and Faure sequences [START_REF] Sobol | On quasi-monte carlo integrations[END_REF]. Sobol low-discrepancy sequences are identified in many aspects as the superior one, compared to other sequences.

Sobol sequences (also called LP τ sequences or (t, s) sequences in base 2) were first introduced by Sobol in 1967 [START_REF] Sobol | On the distribution of points in a cube and the approximate evaluation of integrals[END_REF]. Define a real integrable function f over C n (n-dimensional unit hypercube). The original motivation of Sobol was to generate a sequence x j in C n so that lim

M →∞ 1 M M ∑ j=1 f (x j ) = ∫ C n f (x)dx, (1.64)
and the convergence is as fast as possible.

Sample values of QMC sampling by using Sobol sequences are uniformly placed in input space, even for rather small sample size. Additional conditions, known as property A and A', were also introduced by Sobol in 1975 to make sample values perfectly uniformly distributed in input space [START_REF] Sobol | Uniformly distributed sequences with an additional uniform property[END_REF]:

Definition: A low-discrepancy sequence is said to satisfy Property A if for any binary segment (not an arbitrary subset) of the n-dimensional sequence of length 2 n there is exactly one point in each 2 n hyper-octant that results from subdividing the unit hypercube along each of its length extensions into half.

Definition: A low-discrepancy sequence is said to satisfy Property A' if for any binary segment (not an arbitrary subset) of the n-dimensional sequence of length 4 n there is exactly one point in each 4 n hyper-octant that results from subdividing the unit hypercube along each of its length extensions into four equal parts.

QMC method recently attracts much interest from researchers working in the field of mathematical finance or computational finance. In these fields, high-dimensional numerical integrals in hundreds or thousands of variables are frequently evaluated within a previously defined threshold. Hence, QMC method is widely used in the financial sector, helping value financial derivatives [START_REF] Lecuyer | Quasi-monte carlo methods with applications in finance[END_REF][START_REF] Glasserman | Monte Carlo methods in financial engineering[END_REF].

The detailed generation process of Sobol sequences can be found in Ref. [START_REF] Bratley | Algorithm 659: Implementing sobol's quasirandom sequence generator[END_REF]. An strategies could be found in Refs. [START_REF] Gamboa | Statistical inference for sobol pick-freeze monte carlo method[END_REF][START_REF] Yun | An efficient sampling method for variance-based sensitivity analysis[END_REF].

Determination of analysis results

In the determination of uncertainty and sensitivity analysis results, the most fundamental task is to actually explore the mapping from input vector to output variables:

[x j , y(x j )], j = 1, 2, • • • , M . y

Correlation measures

Correlation indices quantify the correlation and independence between two or more values in fundamental statistics. One of the most often considered correlation indices is Pearson correlation coefficient (CC). The CC was introduced by Pearson in the 1880s, following a related idea of Galton [START_REF] Galton | Opening address to the anthropological section of the british association by the president of the section[END_REF][START_REF] Galton | Regression towards mediocrity in hereditary stature[END_REF][START_REF] Pearson | Note on regression and inheritance in the case of two parents[END_REF]. It provides a measure of the strength of linear correlation between individual input factor and the output variable. The CC between two variables has values between -1 and 1. Three particular values -1, 0, and 1 separately correspond to the total negative linear correlation, no linear correlation, and total linear correlation. Specifically, the CC between an input factor x i and output variable y, denoted by ρ(x i , y), is defined by

ρ(x i , y) = 1 M ∑ M j=1 (x i j -µ i )(y j -E(y)) σ i σ y , (1.65)
where µ i denotes the mathematical expectation of input x i , E(y) the mathematical expectation of output y, σ i the standard deviation of x i , and σ y the standard deviation of y. They are determined by

µ i = 1 M M ∑ j=1 x i j , E(y) = 1 M M ∑ j=1 y j , (1.66)
σ i =   1 M M ∑ j=1 (x i j -µ i ) 2   1/2 , σ y =   1 M M ∑ j=1 (y j -E(y)) 2   1/2
.

(1.67)

The absolute value of ρ(x i , y), between 0 and 1, corresponds to a trend from no linear relationship to an exact linear relationship between x i and y. Typically, for linear (purely additive) models, the sum of squared pairwise CCs between each of the input factors and output variable is 1 in the absence of correlations among input factors: is introduced based on the elements of x except x i :

n ∑ i=1 ρ 2 (x i , y) = 1. ( 1
xi = c 0 + n ∑ j=1,j̸ =i c j x j (1.69)
And the effects of the remaining elements of x on y are represented as

ŷ = b 0 + n ∑ j=1,j̸ =i b j x j .
(1.70)

The PCC between x i and y is defined as the CC between (x i -xi ) and (y -ŷ). An example of the PCCs for nonlinear models will be presented later, together with a similar measure: the partial rank correlation coefficient.

Spearman's rank correlation coefficient (RCC), named after Charles Spearman, characterises the statistical association between the ranking of two variables of interest.

It quantifies how well the dependence between two variables can be explained by a monotonic function. The RCC between x i and y is defined as the CC between their corresponding ranks. Consider a size of sample M . The values of x i and y are replaced by their rank numbers, positive integers between 1 and M , in the sorting of x i and y from the smallest to the largest, namely rx i and ry, respectively. The RCC between x i and y is then given by

ρ(x i , y) = ∑ M j=1 (rx i j -rx i )(ry j -ry) [ ∑ M j=1 (rx i j -rx i ) 2 ] 1/2 [ ∑ M j=1 (ry j -ry) 2 ] 1/2 . (1.71)
When M rank numbers (integers from 1 to M ) are distinct from each other for both x i and y, the mathematical expectation and statistical variance of rank transformed data is calculated as where the formula ∑ M j=1 j 2 = M (2M + 1)(M + 1)/6 was used. The RCC between x i and y then can be simplified as

rx i = ry = M + 1 2 , V (rx i ) = V (ry) = M 2 -1 12 , (1.72)
ρ(x i , y) = 1 - 6 ∑ M j=1 (rx i j -ry j ) 2 M (M 2 -1)
.

(1.73)

The RCC between x i and y has values between -1 and 1, analogous to the CC, with positive values corresponding to the increase or decrease of x i and y along the same direction and negative values corresponding to the increase or decrease of x i and y along opposite directions. A RCC of zero value indicates an obvious non-monotonic relationship between x i and y. The larger the RCC in magnitude, the closer the relationship between x i and y gets to be monotonic. The RCCs are displayed in table 1.3 for three models with scatter plots presented in Figs. 1.10, 1.11 and 1.12, by the use of LHS-based method and one-at-a-time method.

Another rank correlation coefficient is the Kendall rank correlation coefficient(KRCC), commonly referred to as Kendall's τ coefficient. The KRCC was developed by Kendall in 1938 [START_REF] Kendall | A new measure of rank correlation[END_REF]. It measures the ordinal association between two variables under study. The KRCC between x i and y is also defined based on their rank numbers in the sorting from the smallest to the largest. It has values between -1 and 1, with high values indicating a similar (or identical when the coefficient is 1) rank between x i and y and low values labeling a dissimilar (or fully different when the coefficient is -1) rank between them.

A KRCC of zero value between x i and y corresponds to the independent relationship between them. For a size of sample M , an arbitrary pair of mappings from x i to y:

[x i j , y j ], [x i k , y k ], j ̸ = k and j, k = 1, 2, • • • , M
, is said to be concordant if the ranks for both elements follow the same direction, that is, if both x i j > x i k and y j > y k or if both x i j < x i k and y j < y k . It is said to be discordant if the ranks for both elements follow opposite directions, that is, if x i j > x i k and y j < y k or if x i j < x i k and y j > y k . When x i j = x i k or y j = y k , the correspond pair of mappings is neither concordant nor discordant. The definition of KRCC is given by

τ (x i , y) = 2 [(number of concordant pairs)-(number of discordant pairs)] M (M -1)
.

(1.74) In the absence of coupling terms involving different input factors, rank correlation coefficients, including RCC and KRCC, based on one-at-a-time method is equivalent to the PCC and PRCC based on sampling strategies. They all provide reliable measures in quantifying the monotonic relationship between each individual input factor and the output variable. However, in the presence of coupling items, the PRCC is superior to any other correlation coefficient in assessing how well the relationship between an individual input factor and the output variable can be described by a monotonic function.

Occasionally, multiple correlation is also of interest to be considered. The coefficient of multiple correlation, denoted as R in general, assesses how well a given output variable can be predicted using a linear function of a set of input variables [START_REF] Abdi | Multiple correlation coefficient. Encyclopedia of measurement and statistics[END_REF].

R = [ ρ 2 (x 1 , y) + ρ 2 (x 2 , y) -2ρ(x 1 , y)ρ(x 2 , y)ρ(x 1 , x 2 ) 1 -ρ 2 (x 1 , x 2 ) ] 1/2 , (1.75)
where ρ indicates the pairwise CC (Eq. (1.65)) between involved variables. By assuming input variables to be independent of each other, the above expression can be simplified as

R = [ ρ 2 (x 1 , y) + ρ 2 (x 2 , y) ] 1/2 , (1.76)
which is naturally extended to a general situation of n input factors of absent input correlations:

R = [ n ∑ i=1 ρ 2 (x i , y) ] 1/2 . (1.77)
The multiple correlation coefficient is 0.99 for model 1 (Fig. 1.10(a)), 0.84 for model 2 (Fig. 1.11(a)) and 0.89 for SIR model (Fig. 1.12(a))

Sensitivity indices

The sensitivity indices (or measures), also called importance indices (measures), of output variables with respect to input factors are often defined by variance-based methods.

One commonly considered variance-based method is the Sobol's variance decomposition [START_REF] Sobol | Sensitivity estimates for nonlinear mathematical models[END_REF]. Following his concept, the determination of sensitivity indices depends upon multidimensional integrals, as introduced in Section 1.3. Considering the complexity in evaluating multidimensional integrals, the sampling-based method is frequently applied to the calculation of sensitivity indices in practice. With sampling-based method, the sensitivity indices of different orders are evaluated by using two independently constructed samplings, namely A and B, which are M × n matrices. The mathematical expectation and statistical variance of output y are approximated by any of the constructed samplings, for instance, A, as

E(y) = 1 M M ∑ j=1 f (A) j , (1.78)
V (y) = 1 M M ∑ j=1 f 2 (A) j -E 2 (y).
(1.79) Following the concept of Sobol's variance decomposition, the first-order and total sensitivity indices are estimated, respectively, by

     , B (1)
A =      0.
     , A (2)
B =      0.
     , B (2)
A =      0.
     , A (3)
B =      0.
     , B (3)
A =      0.
s i =   1 M M ∑ j=1 f (A) j f (B (i) A ) j -E 2 (y)   / V (y), (1.80)
s T i =   1 M M ∑ j=1 f (A) j ( f (A) j -f (A (i) B ) j )   / V (y), (1.81)
where

A (i)
B is a matrix with column i in A substituted by the ith column in B, and vice versa for matrix B (i)

A whose column i comes from matrix A and the remaining (n -1) columns come from matrix B [START_REF] Saltelli | Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index[END_REF]. An example of construction of both A Analogously, the second-and higher-order sensitivity indices are determined by

s i 1 i 2 =   1 M M ∑ j=1 f (A) j f (B (i 1 i 2 ) A ) j -E 2 (y)   / V (y) -s i 1 -s i 2 , (1.82)
s i 1 i 2 i 3 =   1 M M ∑ j=1 f (A) j f (B (i 1 i 2 i 3 ) A ) j -E 2 (y)   / V (y) -s i 1 -s i 2 (1.83) -s i 3 -s i 1 i 2 -s i 1 i 3 -s i 2 i 3 , (1.84)
. . .

where the first two indices separately quantify the fraction of total variance of output y that is contributed by the interaction effect between x i 1 and x i 2 , and by the interaction effect involving x i 1 , x i 2 and x i 3 .

Alternative formulas for the determination of the first-order and total sensitivity indices are offered by Jansen in 1999 [START_REF] Jansen | Monte carlo estimation of uncertainty contributions from several independent multivariate sources[END_REF][START_REF] Jansen | Analysis of variance designs for model output[END_REF], which are expressed as QMC sampling is much more efficient for sensitivity analysis of models. It provides accurate analysis results by using a small number of sampled points and also avoids highly computational cost in performing configuration average.

s i = 1 -   1 2M M ∑ j=1 ( f (B) j -f (A (i) B ) j ) 2   / V (y), (1.85)
s T i =   1 2M M ∑ j=1 ( f (A) j -f (A (i) B ) j ) 2   / V (
               , A =                0.
               . B =                0.
               . Chapter 2
The analytic analysis for models with independent input variables

Any operations that we perform on a model response dependent upon a number of input variables of uncertainty require us to identify the response uncertainty based on the uncertainty in input variables. The propagation of variance, characterising the effect of input uncertainty on the uncertainty of model response, constitutes the essential ingredient of uncertainty and sensitivity analysis of complex systems. In the present chapter, an analytic formula is derived by using Taylor series to specify the variance propagation from input variables to the model response in the absence of input correlations. With the formula, we can exactly calculate the uncertainty of model response. Furthermore, it also allows one to quantify the partial variance contributions of different orders from input variables to the output one, whereby input variables can be ranked according to their importance in explaining the uncertainty of the output variable.

Taylor series

The concept of a Taylor series was formulated by a Scottish mathematician James Gregory and formally introduced by an English mathematician Brook Taylor in 1715. In mathematics, a Taylor series is a representation of a function as an infinite sum of terms that are evaluated from the values of the function's derivatives at a single point [START_REF] Taylor | Direct and Reverse Methods of Incrementation" Translated into English in[END_REF].

In general, a function can be approximated as the sum of a finite number of terms constituting of function's derivatives of different orders by using Taylor series of neglecting higher-order terms.

Taylor series is widely used for the expansion of functions only dependent upon one single variable. Regarding a general univariate function of the form y = f (x), its Taylor series at a point x = a is given by

∞ ∑ i=0 f (i) (a) i! (x -a) i , (2.1)
where i! denotes the factorial of i, and f (i) (a) the ith derivative of f evaluated at point a.

The zeroth-order derivative of f is defined as f itself and 0!=1. If a = 0, the expansion is known as a Maclaurin series, named after the Scottish mathematician Colin Maclaurin who introduced extensive use of this special case of Taylor series in the 18th century [START_REF] Grabiner | Was newton's calculus a dead end? the continental influence of maclaurin's treatise of fluxions[END_REF].

The Taylor series can also be generalised to functions of more than one input variable [START_REF] Hörmander | The analysis of linear partial differential operators I: Distribution theory and Fourier analysis[END_REF]. For a real function of the form

y = f (x) with x = (x 1 , x 2 , • • • , x n ) T labeling
the input vector of n-dimensional variables, its Taylor series at a specified point

x = (a 1 , a 2 , • • • , a n ) T is represented as ∞ ∑ i 1 =0 • • • ∞ ∑ in=0 (x 1 -a 1 ) i 1 • • • (x n -a n ) in i 1 ! • • • i n ! ( ∂ i 1 +•••+in f ∂x i 1 1 • • • ∂x in n ) (a 1 , • • • , a n ). (2.2)

Variance propagation for univariate case

Beginning with the univariate function y = f (x), it can be expanded by the use of Taylor series at the central point of x as

y = f (µ) + ∞ ∑ i=1 1 i! ( d i f dx i x=µ ) (x -µ) i , (2.3)
in which µ indicates the mathematical expectation of x. Taking the average of both sides of Eq.(2.3) yields

E(y) = f (µ) + ∞ ∑ i=1 1 i! ( d i f dx i x=µ ) M i (x), (2.4)
where M i (x) is the i th central moment of variable x, given by

M i (x) = ∫ (x -µ) i P (x)dx.
(2.5) P (x), as stated in Chapter 1, denotes the PDF of x and characterises the epistemic uncertainty in x. The variance of y, namely V (y), is defined in mathematics by

V (y) = ∫ y 2 P (x)dx -E 2 (y). (2.6)
Substituting Eqs. (2.3) and (2.4) into Eq. (2.6) provides

V (y) = ∫ [ f (µ) + ∞ ∑ i=1 1 i! ( d i f dx i x=µ ) (x -µ) i ] 2 P (x)dx - [ f (µ) + ∞ ∑ i=1 1 i! ( d i f dx i x=µ ) M i (x) ] 2 =f 2 (µ) + ∞ ∑ i,j=1 1 i! • j! ( d i f dx i x=µ • d j f dx j x=µ ) ∫ (x -µ) i+j P (x)dx + 2f (µ) ∞ ∑ i=1 1 i! ( d i f dx i x=µ ) ∫ (x -µ) i P (x)dx -2f (µ) ∞ ∑ i=1 1 i! ( d i f dx i x=µ ) M i (x) - ∞ ∑ i,j=1 1 i! • j! ( d i f dx i x=µ • d j f dx j x=µ ) M i (x)M j (x) -f 2 (µ) = ∞ ∑ i,j=1 1 i! • j! ( d i f dx i x=µ • d j f dx j x=µ ) [M i+j (x) -M i (x)M j (x)] , (2.7)
which is equivalent to

V (y) = ∞ ∑ i,j=0 1 i! • j! ( d i f dx i x=µ • d j f dx j x=µ ) [M i+j (x) -M i (x)M j (x)] , (2.8)
since M 0 (x) = 1. This analytic expression mathematically explains the output variance propagated from the uncertainty in input variable through the functional relationship connecting y and x. Furthermore, it can also help identify the partial variance contributions of different orders of x, embodied by the derivatives of different orders of f . If one just considers the first-order variance contribution of x, Eq. (2.8) is then simplified as

V (y) ≈ ( df dx x=µ ) 2 V (x), (2.9)
where M 1 = 0 and M 2 = V (x) have been used, and V (x) labels the variance of x.

Equation (2.9), called the general Taylor series expansion truncated to the first order, is most widely used for the estimation of the uncertainty in y in terms of the uncertainty in x. This approximation, however, is satisfactory for the uncertainty analysis of highly nonlinear models only when the uncertainty in input variable is negligible [START_REF] Padulo | Comparative analysis of uncertainty propagation methods for robust engineering design[END_REF], as stated in the first chapter.

To identify the partial variance contributions of different orders from input x, a new quantity of interest is proposed, labeled as R γ v . It is defined as the ratio of V (y) contributed by first γ orders of uncertainty in x to its exact value:

R γ v = V γ (y) V T (y) , (2.10)
where V γ (y) is calculated with Eq. (2.8) under the condition i + j ≤ 2γ and V T (y) the theoretical value of V (y) obtained from the below integral:

V T (y) = ∫ y 2 P (x)dx -E 2 (y), = ∫ f 2 (x)P (x)dx - (∫ f (x)P (x)dx ) 2 .
(2.11)

The analytic formula previously deduced (Eq. (2.8)) specifies the variance propagation from input variable to the model output. With this formula, it is natural to quantify the importance of partial variance contributions of different orders from input x in the estimation of uncertainty in output y. In the following, higher-order contributions of input uncertainty to the uncertainty of output y are evaluated for highly nonlinear models, by using the analytic formula. Two widely discussed distribution laws: uniform distribution and normal distribution, as examples, specify the uncertainty in input x.

Uniform distribution

We first suppose x is uniformly distributed. Its PDF is given by

P (x) =        0, x < x 0 , 1 xm-x 0 , x 0 ≤ x ≤ x m , 0, x > x m .
(2.12)

Uniform distribution provides the mathematical expectation of x: µ = xm+x 0 2 and its central moments:

M k (x) = { 0, k is odd, 3 k/2 k+1 σ k , k is even, (2.13)
where k is a positive integer; σ labels the standard deviation of x; yields

V (y) = ∞ ∑ i=1 ∞ ∑ j=i j=j+2 C ij ( d i f dx i x=µ • d j f dx j x=µ ) σ i+j , (2.14)
where j is summed with an increment of 2 and C ij is defined as follows.

If i = j C ij =      3 (i+j)/2
(i+j+1)•i!•j! , i and j are odd,

3 (i+j)/2 •i•j (i+j+1)•(i+1)!•(j+1)!
, i and j are even.

(2.15)

Else

C ij =      2•3 (i+j)/2 (i+j+1)•i!•j! , i and j are odd, 2•3 (i+j)/2 •i•j (i+j+1)•(i+1)!•(j+1)!
, i and j are even.

(2.16)

The underlying results of R γ v are displayed in Fig. 2.1, for two widely discussed nonlinear functions: power-law function and exponential function. It is noticed that R γ v will be close to one while γ is large enough for both kinds of functions with different parameters, see left plots in both panels.

By inserting the functions under discussion: y = x α and y = e αx , equation (2.14) can be separately updated by

∞ ∑ i=1 ∞ ∑ j=i j=j+2 C ij (α!) 2 (α -i)! • (α -j)! µ 2α ( σ µ ) i+j , (2.17)
and

∞ ∑ i=1 ∞ ∑ j=i j=j+2 e 2αµ C ij (ασ) i+j , (2.18)
which demonstrate that, for specified γ (thereby both i and j are specified) and α, the contribution of uncertainty in input x to the observation R γ v is embodied in the ratio of σ to µ for power-law relationships connecting y and x and in σ for exponential relationships. In Fig. 2.1, middle plots of both panels display the distribution of R γ v versus exponent α under specified µ/σ (for power-law function) and fixed σ (for exponential function), by considering different order contributions of uncertainty in x. Right ones illustrate the dependence of R γ v on µ/σ (1/σ) for power-law function (exponential function) by, particularly, setting α = 4. Regarding power-law function, γ should be large for large |α -1| if µ/σ is constant, or for small µ/σ if α is constant, in order to make R γ v more closer to 1, i.e., to make the estimated variance in model output more closer to its exact value. For example, the contribution of uncertainty in x should be considered until up to the 6th order (to ensure R γ v > 0.98) in regions of α < 0.1 and α > 8. For the exponential relationship connecting y and x, R γ v is symmetric with respect to α = 0. The contributions of higher orders of input uncertainty are non-ignorable when the functional form deviates from linear law and also σ is non-negligible. The statement is visually verified that the original approximation presented in Eq. (2.9), by only considering the contribution of the first order of input uncertainty, can successfully estimate the uncertainty in output variable only when the input uncertainty is negligible or the model under discussion is almost linear. The contributions of higher order(γ ≥ 2) of input uncertainty, however, are non-ignorable for models behaving non-linearly in the neighbourhood of input variable when σ (for the exponential association connecting input and output variables) or σ/µ (for the power-law association connecting input and output variables) is non-negligible.

Normal distribution

Optionally, specify the uncertainty in x by normal (Gaussian) distribution. Its PDF is given by which provides

P (x|µ, σ) = 1 σ √ 2π e -(x-µ) 2 2σ 2 ; x ∈ (-∞, +∞), (2.19)
M k (x) = { 0, k is odd, σ k (k -1)!!, k is even, (2.20)
where (⋆)!! is the double factorial with definition given by

(⋆)!! =        ⋆ • (⋆ -2) • (⋆ -4) • • • 5 • 3 • 1, ⋆ > 0 odd, ⋆ • (⋆ -2) • (⋆ -4) • • • 6 • 4 • 2, ⋆ > 0 even, 1, ⋆ = -1, 0.
(2.21)

Associations of the observation R γ v with order index γ, exponent parameter α, as well as the input uncertainty that is specified by µ and σ, are presented in Fig. 2 displays the results for the exponential relationship. The observation R γ v will be close to one when γ is large enough for each constant α, see left plots in both panels. We can notice from the middle plots that the contributions of higher orders of input uncertainty should be considered for highly nonlinear functions. Right plots in both panels suggest that R γ v of γ = 1 can reach 1 only when the ratio of µ to σ is very large for power-law function (y = x 5 as an example) or σ is very small for exponential function (y = e x as an example). Typically, when the input variable is normally distributed, parameter α of power-law function should be a positive integer to enable the mean and variance of output y calculable.

To understand the influence of different distribution laws specifying the uncertainty in input variable on the estimation of output variance, a comparison of uniform and normal distributions is provided for the observation R γ v by keeping both µ and σ the same for two kinds of law, see Fig. 2.3. The exponential relationship connecting y and x: y = e αx , as an example, is discussed here. Results suggest that, when both α and σ are specified, uniformly distributed input variable will drive the evaluated output uncertainty more close to its exact value, compared to normally distributed one, in considering variance contributions of the same order.

Generalisation of the analytic formula

The analytic formula Eq. (2.8) for variance propagation is valid only for the model of a single input variable. However, in mathematical modeling, the output quantity often depends upon two or more input variables of uncertainty. In the following, the analytic formula is generalised to the general case of more than one input variable. Recall the real model of n independent input variables of the form

y = f (x) = f (x 1 , x 2 , • • • , x n ). (2.22)
By the use of multi-variate Taylor series, the model response y can be expanded at the

central point ({µ}) = (µ 1 , µ 2 , • • • , µ n ) of input vector as y = f ({µ}) + n ∑ t=1 ∞ ∑ i=1 (x t -µ t ) i i! ( ∂ i f ∂x i t )({µ}) + n ∑ t,s=1 t<s ∞ ∑ it,is=1 (x t -µ t ) it (x s -µ s ) is i t ! • i s ! × ( ∂ it+is f ∂x it t ∂x is s )({µ}) + • • • + ∞ ∑ i 1 ,••• ,in=1 (x 1 -µ 1 ) i 1 • • • (x n -µ n ) in i 1 ! • • • i n ! ( ∂ i 1 +•••+in f ∂x i 1 1 • • • ∂x in n )({µ}), (2.23)
where µ i denotes the mathematical expectation of input x i , as before. Reviewing the mathematical definition of variance, V (y) is analogously calculated as

V (y) = ∞ ∑ i 1 •••in=0 j 1 •••jn=0 F (x 1 ) i 1 j 1 •••(xn) injn A(i 1 , • • • , i n , j 1 , • • • , j n ) ( ∂ i 1 +•••+in f ∂x i 1 1 • • • ∂x in n • ∂ j 1 +•••+jn f ∂x j 1 1 • • • ∂x jn n ) ({µ}), (2.24)
with

A(• • • ) = i 1 ! • • • i n ! • j 1 ! • • • j n !, (2.25)
and

F (x 1 ) i 1 j 1 •••(xn) injn = M i 1 +j 1 (x 1 ) • • • M in+jn (x n ) -M i 1 (x 1 )M j 1 (x 1 ) • • • M in (x n )M jn (x n ).
(2.26) Equation (2.24) specifies the variance propagation from input variables to the output one for a general model, by assuming input variables are statistically independent of each other.

Applications to the analysis of complex systems

The power grid system and economic system are considered in this section as examples of the application of our generalised analytic formula, Eq. (2.24), in the uncertainty and sensitivity analysis of complex systems. These two systems play extremely important roles in modern society. Their reliability analyses have attracted many researchers' interest.

Concerning the topic of sensitivity analysis, someone is of most interest to the sensitivity indices. They allow one to quantitatively interpret the importance of individual factors and of their interaction effects in the estimation of model response. Recalling the analytic formula of variance propagation in general models (see Eq. (2.24)), partial variance contributions of different dimensionality in Eq. (1.33) can be computed through

V t = ∞ ∑ i,j=0 F (xt) ij A(i, j) ( ∂ i f ∂x i t × ∂ j f ∂x j t )({µ}), (2.27)
V ts = ∞ ∑ i,j,k,l=0 F (xt) ij (xs) kl A(i, j, k, l) ( ∂ i+k f ∂x i t ∂x k s • ∂ j+l f ∂x j t ∂x l s ) ({µ}) -V t -V s , (2.28)
V tsu = ∞ ∑ i,j,k,l,p,q=0 F (xt) ij (xs) kl (xu) pq A(i, j, k, l, p, q) ( ∂ i+k+p f ∂x i t ∂x k s ∂x p u • ∂ j+l+q f ∂x j t ∂x l s ∂x q u ) ({µ}) -V t -V s -V u -V ts -V tu -V su , (2.29)
. . .

The sensitivity indices of different orders are then quantified accordingly by

s t = V t V (y) , s ts = V ts V (y) , s tsu = V tsu V (y) , • • • , s T t = s t + n ∑ s=1;s̸ =t s ts + n ∑ s,k=1;s>k;s,k̸ =t s tsk + • • • . (2.30)

Power grid system

In our daily life, the main form of energy that we think of is the power from electricity.

Typically, a complex distribution system used for the transmission of electric power is called the power grid. The electricity power grid is a physical system that delivers electricity from the place where it is generated to the site where it is used. Consequently, a power grid system can be segmented into three sections: electric power generation, supply, and transmission, which work together to meet the basic electricity demand of ordinary people. Uncertainty and reliability analysis of the power grid system has been carried out since 1994 by using MC methods [START_REF] Billinton | Reliability Assessment of Electric Power Systems Using Monte Carlo Methods[END_REF]. In general, the assessment of power grid system reliability is divided into two aspects: system adequacy and system security which are, respectively, related to the steady-state operation of system and to the ability of system to withstand sudden natural disturbances or to avoid attack.

Facing the increased global energy consumption, electric utility companies are striving to generate wind power to meet the growing electricity demand. Wind power, as a typical example of renewable energy sources, is plentiful, widely distributed, clean, produces no greenhouse gas emissions during operation, consumes no water, and uses little land [START_REF] Fthenakis | Land use and electricity generation: A life-cycle analysis[END_REF]. In this part, a framework on the uncertainty and sensitivity analysis of the actual wind power output, denoted as P d , is established by using the above established analytic method.

In theory, the wind output power of the wind turbine generator is always the same when the wind speed is specified. In reality, however, the output power for a fleet wind turbine generator of the same type exhibits considerable variations even when they are operating at the same wind speed [START_REF] Thiringer | Control by variable rotor speed of a fixed-pitch wind turbine operating in a wide speed range[END_REF][START_REF] Slootweg | Modeling of large wind farms in power system simulations[END_REF]. Consequently, P d is often considered as a random variable that is characterised by two parameters: x and ε through a functional relationship:

P d (x) = f (x) + ε, (2.31)
where x indicates the wind speed, f (x) the deterministic power output from a wind turbine generator, which acts as a function of the wind speed:

f (x) =                  0, 0 ≤ x ≤ v ci , (A + Bx + cx 2 ) * P r , v ci ≤ x ≤ v r , P r , v r ≤ x ≤ v co , 0, x ≥ v co , (2.32)
and ε the variation of the power output, obeying normal distribution with µ = 0 and σ = 0.1P r [START_REF] Park | Planning manual for utility application of wecs[END_REF][START_REF] Giorsetto | Development of a new procedure for reliability modeling of wind turbine generators[END_REF][START_REF] Jin | Uncertainty analysis for wind energy production with dynamic power curves[END_REF][START_REF] Sansavini | A stochastic framework for uncertainty analysis in electric power transmission systems with wind generation[END_REF]. v ci is called the cut-in wind speed at which the turbine first starts to rotate and generate power, v r the rated wind speed at which the rated power P r (the power output limit that the electrical generator is capable of) is reached, and v co the cut-out wind speed at which the wind turbine generator is shut down for safety reasons, as shown in Fig. 2.4. Following the reference [START_REF] Sansavini | A stochastic framework for uncertainty analysis in electric power transmission systems with wind generation[END_REF], v ci = 3ms -1 and v r = 12ms -1 . The constants A, B, and C are determined by v ci and v r through equations [219] In the power grid system, the uncertainty of the wind speed x is often represented by a

A = 1 (v ci -v r ) 2 [ v ci (v ci + v r ) -4v ci v r ( v ci + v r 2v r ) 3 ] , (2.33)
B = 1 (v ci -v r ) 2 [ 4(v ci + v r ) ( v ci + v r 2v r ) 3 -(3v ci + v r ) ] , (2.34)
C = 1 (v ci -v r ) 2 [ 2 -4 ( v ci + v r 2v r ) 3 ] . ( 2
Weibull distribution of two parameters as

P (x) = k c ( x c ) k-1 e -( x c ) k ; (k > 0, x > 0, c > 1), (2.37)
where c and k are the scale parameter and the shape parameter, respectively, [START_REF] Johnson | Wind energy systems[END_REF].

Three different combinations of c and k are considered here. The underlying probability density function of the wind speed is displayed in Fig. 2.5. In our analysis, the wind speed x is limited to the real range [3m/s, 12m/s]. Parameters characterising the uncertainty in x are displayed in table 2.1. The exact values obtained by our method are in good agreement with those provided by Sobol's method. This implies our method is equivalent to Sobol's one in the analysis of the model under consideration.

Economic system

Regarding the economic system, one of the oldest classic production scheduling models is the EOQ model. This model was developed by Harris in 1913 [START_REF] Harris | How many parts to make at once[END_REF]. It aims at determining the economic order quantity that minimizes the total system cost. Some analyses about the uncertainty and sensitivity of this model have been proposed in Refs. [START_REF] Borgonovo | Sensitivity analysis: A review of recent advances[END_REF][START_REF] Borgonovo | Global sensitivity analysis in inventory management[END_REF][START_REF] Schwarz | The economic order-quantity (eoq) model[END_REF]. However, the discussion of its reliability to each input parameter, especially to the interaction effects between different parameters, is still limited.

This subsection builds an intuitive insight into the uncertainty and reliability of EOQ model in terms of the analytic formula deduced before. In EOQ model, three types of cost are accounted: purchase cost, holding cost and ordering cost. The purchase cost is the variable cost of goods, equaling the product of the price of the units themselves, denoted as P , and the annual demand quantity, denoted as D. The price of the units themselves is assumed to be fixed, regardless of the number of demand quantity. The holding cost indicates the cost of holding units in inventory (management's cost of capital invested in the units, the cost of the space consumed by the units, taxes of insurance, etc). It is expressed as the product of the cost of holding unit, denoted as h, and the average quantity in inventory (between fully replenished and empty), denoted as Q/2. Q

is the order quantity. The ordering cost represents all the costs associated with placing orders excluding the purchase cost. Suppose each order has a fixed cost K. The number of orders we need to make per year is D/Q. Hence the ordering cost is K * D/Q. The total system cost in EOQ model can then be represented as

T C = P D + DK Q + hQ 2 . (2.38)
EOQ is the order quantity that minimizes the total system cost. It is easy to obtain the value of Q which determines the minimum point of T C, that is

Q * = √ 2DK h . (2.39)
The uncertainty of Q * , as well as its sensitivity with respect to independent input variables D, K, h and to their interactions, is quantified in this part by using our analytic method. the analysis results are almost equal to the exact values. This suggests that the contributions of input uncertainty of third-or higher-order can be neglected and the analysis results of γ = 2 can truly represent the reliability of EOQ model. Results state that all three parameters are important in establishing the uncertainty in model response.

Furthermore, the interaction effect between each two parameters also contributes a small part to the uncertainty in response Q * . s DK = -0.002 says that the interaction effect between D and K will result in a small decrease in the variance of model response.

A comparison of our analysis results with Sobol's values is shown in Fig. 2.7. The convergent results of our method are almost the same as Sobol's numbers.

Chapter 3

The analytic analysis for models with correlated input variables

Over the past few decades, the uncertainty and sensitivity analysis of models is mainly performed by assuming that input variables are independent of each other. A problem is often arisen in the development of the methodology for sensitivity analysis: how to interpret the sensitivity measures when input variables are non-independent? Recently, the interest in uncertainty and sensitivity analysis is increasing in the presence of correlated input variables, facing the existence of correlated factors in practical applications.

In this chapter, we focus on the extension of the analytic framework established in the previous chapter. The extended framework allows one to understand, from a theoretical point of view, the variance propagation from correlated input variables to the model response and also the sensitivity of model response with respect to input independence and correlations.

Variance propagation

Analogously, consider a general model of the form y

= f (x) with x = (x 1 , x 2 , • • • , x n ) T
labeling the input vector of n-dimensional variables. Recalling the multi-variate Taylor series at the central point of input vector, presented in Eq. (2.23), the mathematical expectation of the model response can be represented in the presence of input correlations as

E(y) =f ({µ}) + n ∑ t=1 ∞ ∑ it=1 M it (x t ) A(i t ) ( ∂ it f ∂x it t )({µ}) + n ∑ t,s=1 t<s ∞ ∑ it,is=1 M it,is (x it t , x is s ) A(i t , i s ) ( ∂ it+is f ∂x it t ∂x is s )({µ}) + • • • + ∞ ∑ i 1 ,••• ,in=1 M i 1 ,••• ,in (x i 1 1 , • • • , x in n ) A(i 1 , • • • , i n ) ( ∂ i 1 +•••+in f ∂x i 1 1 • • • ∂x in n )({µ}) = ∞ ∑ i 1 ,••• ,in=0 M i 1 ,••• ,in (x i 1 1 , • • • , x in n ) A(i 1 , • • • , i n ) ( ∂ i 1 +•••+in f ∂x i 1 1 • • • ∂x in n )({µ}), (3.1)
where A(

M i 1 ,••• ,in (x i 1 1 , • • • , x in n ) = E [ (x 1 -µ 1 ) i 1 • • • (x n -µ n ) in ] = ∫ (x 1 -µ 1 ) i 1 • • • (x n -µ n ) in P (x)dx, (3.2)
with P (x) indicating the joint PDF of x. In the absence of input correlations, P (x) is simplified as

P (x) = n ∏ i=1 P i , (3.3)
where P i denotes the PDF of input x i , as before. In the existence of correlations among input variables, the variance of model response is derived by the use of Eqs. (2.23) and

(3.1) as

V (y) = ∞ ∑ i 1 ,••• ,in=0 j 1 ,••• ,jn=0 1 A(i 1 , • • • , i n , j 1 , • • • , j n ) ( ∂ i 1 +•••+in f ∂x i 1 1 • • • ∂x in n • ∂ j 1 +•••+jn f ∂x j 1 1 • • • ∂x jn n ) ({µ}) × [ M i 1 +j 1 ,••• ,in+jn (x i 1 +j 1 1 , • • • , x in+jn n ) -M i 1 ,••• ,in (x i 1 1 , • • • , x in n )M j 1 ,••• ,jn (x j 1 1 , • • • , x jn n )
] .

(3.4)

Apparently, in the absence of input correlations, the above equation is coincident with Eq. (2.24) that explains the variance propagation from independent input variables to the model response.

The concept of complete variance decomposition, presented in Eq. (1.33), is proposed by assuming input independence. Its form is also valid for the correlated case. In the presence of correlated input variables, however, partial variance contributions with dimensionality larger than 1 are contributed not only by the coupling items presented in the functional form of the model under discussion (for independent case), but also by the input correlations. Regarding the situation of correlated inputs, the impact of a single variable can be represented as the sum of contributions provided by its correlations with the remaining variables and by its independence, respectively. Based on the description, each fractional variance contribution included in the original variance decomposition can be divided into three sections: independent variance contribution (labeled by superscript U), correlated variance contribution (labeled by superscript C), and coupling variance contribution (labeled by superscript UC). Mathematically, the output variance can be decomposed in the presence of input correlations as

V (y) = n ∑ i=1 (V U i + V C i + V UC i ) + n ∑ i,j=1 i<j (V Up ij + V Cp ij + V UCp ij ) + • • • + (V Uq 12•••n + V Cq 12•••n + V UCq 12•••n ), (3.5)
where p ∈ {i, j}, q ∈ {1, 2, • • • , n}, and

V i = V U i + V C i + V UC i , (3.6)
V ij = V Up ij + V Cp ij + V UCp ij , (3.7)
. . .

V 12•••n = V Uq 12•••n + V Cq 12•••n + V UCq 12•••n . (3.8)
V U i (V C i )
is the variance contribution produced by the independent (correlated) section of x i alone, V U i ij (V C i ij ) the contribution of the interaction between x j and the independent (correlated) section of x i , and so on up to

V Uq 12•••n (V Cq 12•••n
) the contribution of the interaction associated with the independent (correlated) section of x q and the rest variables. Coupling variance contributions are produced by the coupling effects between independent and correlated sections for individual input variables.

Estimation of sensitivity indices

Working within a probabilistic framework, variance-based sensitivity measures are defined on the bases of partial contributions presented in the variance decomposition of model response. In the determination of each partial variance contribution, higher-order covariance embodied in the analytic formula Eq. (3.4) should be concerned for nonlinear models. Consequently, it is necessary to specify the correlated and independent sections of single input variables for the confirmation of fractions contained in Eq. (3.5), whereby the importance of the independent section, correlated section and their coupling effect can be quantified for each individual input variable, in establishing the uncertainty of model response.

Generation of correlated variables

In probabilistic models the dependency between input variables is often represented by the CC that is defined as

ρ(x i , x j ) = E[(x i -µ i )(x j -µ j )] σ i σ j , (3.9)
where E[ * ] is the expectation operation by returning the average value of * . Eq. (3.9) is equivalent to Eq. (1.65) which defines the CC between an arbitrary input variable and the model response. For the sake of simplicity in writing, ρ(x i , x j ) is simplified as ρ ij in the following discussion.

In the presence of correlations, an arbitrary variable can be represented as the sum of a correlated section and an independent section. Based on the description, we can rewrite the input vector as

x = x C + x U (3.10) = (x C 1 , x C 2 , • • • , x C n ) T + (x U 1 , x U 2 , • • • , x U n ) T . (3.11)
The correlated sections x C i with i = 1, 2, • • • , n indicate the correlations of x i with the remaining variables. By using the linear correlation model (Eq. (3.9)), x C can be generated as

x C = Ax, (3.12)
where A is an n × n nonsymmetric hollow matrix (diagonal elements are all equal to zero). The independent sections x U i with i = 1, 2, • • • , n denote the independence of x i . They are often represented by newly introduced random variables. x U can be expressed as

x U = Cr, (3.13)
where C is an n × n diagonal matrix and r = (r 1 , r 2 , • • • , r n ) T an n-dimensional vector of random variables. r must follow the same type distributions as x and also satisfy

µ(r) = C -1 (µ -Aµ), (3.14
a ij = σ i σ j     1 - ∑ k<l, k,l̸ =i ρ 2 kl (1 - ∑ h<q, h>k;h,q̸ =i,l ρ 2 hq ) + 2 ∑ k<l<h, k,l,h̸ =i ρ kl ρ kh (ρ lh - ∑ q̸ =i,l,h;q>k ρ lq ρ hq )     -1 ×     ρ ij (1 - ∑ k<l, k,l̸ =i,j ρ 2 kl + 2 ∑ k<l<h, k,l,h̸ =i,j ρ kl ρ kh ρ lh ) - ∑ k̸ =i,j ρ ik ρ jk (1 - ∑ h<q, h,q̸ =i,j,k ρ 2 hq ) + ∑ k<l, k,l̸ =i,j (ρ ik ρ jl + ρ il ρ jk )(ρ kl - ∑ h̸ =i,j,k,l ρ kh ρ lh )     , (3.16)
and

c ii =     1 - ∑ k<l, k,l̸ =i ρ 2 kl (1 - ∑ h<q, h>k;h,q̸ =i,l ρ 2 hq ) + 2 ∑ k<l<h, k,l,h̸ =i ρ kl ρ kh (ρ lh - ∑ q̸ =i,l,h;q>k ρ lq ρ hq )     -1/2 ×     1 - ∑ k<l ρ 2 kl (1 - ∑ h<q, h>k;h,q̸ =l ρ 2 hq ) + 2 ∑ k<l<h ρ kl ρ kh (ρ lh - ∑ q̸ =l,h;q>k ρ lq ρ hq )     1 2
, (3.17)

where the fifth-and higher-order terms are neglected. The above expressions, derived according to the analysis of simple cases as shown in Appendix B, constitute essential ingredients of quantifying sensitivity measures associated with correlated section, independent section, and their coupling effect, for each single input variable.

If two variables, say x i and x j , simultaneously correlate with a third variable, say x k , the correlation between x i and x j is then not arbitrary assigned. With given correlation coefficients ρ ik and ρ jk , x i and x j then can be formulated based on x k as

x i = ρ ik σ i σ k x k + √ 1 -ρ 2 ik u i , (3.18)
and

x j = ρ jk σ j σ k x k + √ 1 -ρ 2 jk u j , (3.19)
where u i and u j are newly introduced variables that are independent of x k . Inserting the above formulations into Eq. (3.9) yields

ρ ij = ρ ik ρ jk + √ (1 -ρ 2 ik )(1 -ρ 2 jk )ρ(u i , u j ). (3.20)
In statistics, ρ(u i , u j ), labeling the linear correlation between variables u i and u j , can be as large as 1 (total positive correlation) and as small as -1 (total negative correlation).

This implies the correlation ρ ij between x i and x j will be limited to the range

[ ρ ik ρ jk - √ (1 -ρ 2 ik )(1 -ρ 2 jk ), ρ ik ρ jk + √ (1 -ρ 2 ik )(1 -ρ 2 jk ) ] (3.21)
when the correlations of them with a third variable x k are specified. The relationship ρ ij = ρ ik ρ jk holds iff the newly introduced variables u i and u j are set to be independent of each other.

Sensitivity indices

With help of the analytic formula (3.4) that explains the variance propagation in the presence of input correlations, the partial variance contributions of different dimensionality can be calculated by

V i = ∞ ∑ k,l=0 1 k! • l! ( ∂ k f x k i • ∂ l f ∂x l i )({µ}) • [M k+l (x i ) -M k (x i )M l (x i )] , (3.22)
V ij = ∞ ∑ k,l,p,q=0 1 k! • l! • p! • q! ( ∂ k+p f x k i x p j • ∂ l+q f ∂x l i x q j )({µ}) [ M k+l,p+q (x k+l i , x p+q j )- M k,p (x k i , x p j )M l,q (x l i , x q j ) ] -V i -V j , (3.23)
. . .

where multivariate central moments with different variables are involved. The (k, l) central moment with any two different variables, say x i and x j , is defined as (simplified from Eq. (3.2))

M k,l (x k i , x l j ) = E [ (x i -µ i ) k (x j -µ j ) l ] . (3.24)
It can be derived analytically by formulating one variable on the basis of another:

x i = ρ ij σ i σ j x j + √ 1 -ρ 2 ij u i , (3.25)
or

x j = ρ ij σ j σ i x i + √ 1 -ρ 2 ij u j , (3.26)
with u i (u j ) independent of x j (x i ) and taking the same standard deviation as x i (x j ). If k ̸ = l, the above two formulating strategies are equivalent in determining M k,l (x k i , x l j ) only when x i and x j are normally distributed and hold the same standard deviation.

The kth-order central moment of x is given by

M k 1 ,••• ,km (x) = E [ m ∏ i=1 (x i -µ i ) k i ] , (3.27)
where

k 1 + k 2 + • • • + k m = k.
For normal distribution, we have

M k 1 ,••• ,km (x -µ) = { 0, k is odd, ∑ (σ ij σ lh • • • ), k is even, (3.28)
where the sum is taken over all allocations of the set {1, 2, • • • , k} into k/2 pairs, and σ ij the covariance of x i and x j . For example, for the fourth-order central moment (k = 4) with four different variables (calculated in detail in Appendix D), one sums the products of any two covariances:

E [ 4 ∏ i=1 (x i -µ i ) ] = σ 12 σ 34 + σ 13 σ 24 + σ 14 σ 23 . (3.29)
This yields (k -1)!! (double factorial) terms in the sum. When nonlinear terms of a variable are contained in the central moments, like (x i -µ i ) k in Eq. (3.24), one can expand it as the product of k variables and uses σ ii = σ 2 i . Consequently, for normally distributed x i and x j with σ i = σ j , M k,l (x k i , x l j ) is invariant in exchanging x i and x j , that is

M k,l (x k i , x l j ) = M l,k (x l i , x k j ). (3.30)
This property (invariant in exchanging any two variables) could go for any order central moments with any number of different variables that follow a normal distribution and have the same standard deviation.

For uniformly distributed variables, however, the fourth-order central moment with four different variables is calculated as

E [ 4 ∏ i=1 (x i -µ i ) ] = σ 12 σ 34 + σ 13 σ 24 + σ 14 σ 23 - 6 5 ∆, (3.31)
where ∆ is a product of three covariances, dependent upon how to choose a variable to be unformulated in the calculating process. For example, when x 2 is selected to be unformulated, we will get ∆ = σ 12 σ 23 σ 24 .

(3.32)

Furthermore, if one focuses on the correlated, independent and coupling effects contained in a high-order covariance, see Eq. (3.24), x i (x j ) should be formulated on the bases of all the remaining input parameters, as discussed before.

The total contributions to the variance of model response, provided by the independence, correlation, and their coupling effect are represented, for an arbitrary variable x i , as

V TU i = V U i + n ∑ j̸ =i V U i ij + • • • + V U i 12•••n , (3.33)
V TC i = V C i + n ∑ j̸ =i V C i ij + • • • + V C i 12•••n , (3.34)
V TUC i = V UC i + n ∑ j̸ =i V UC i ij + • • • + V UC i 12•••n . (3.35)
The sensitivity (or importance) measures can then be determined by

s U i = V U i V (y) , s C i = V C i V (y) , s UC i = V UC i V (y) , s TU i = V TU i V (y) , s TC i = V TC i V (y) , s TUC i = V TUC i V (y) . (3.36)
The first three measures are called the main sensitivity indices which, separately, denote the importance of the independent section, correlated section, and their coupling effect, by neglecting the interaction effects among input parameters. The last three measures are called the total sensitivity indices which denote the importance of each corresponding part in considering the interaction effects of x i with the remaining inputs. Similar to Sobol's indices, the main and total sensitivity indices are defined as

s i = V i V (y) = s U i + s C i + s UC i , (3.37)
S T i =   V i + ∑ j̸ =i V ij + • • •   / V (y) = s TU i + s TC i + s TUC i , (3.38)
which evaluate the importance of individual input parameters x i before and after considering the interaction effects among different input parameters, respectively.

With given correlation coefficients between any two input parameters, our method allows one to evaluate the importance of individual input parameters in the estimation of model output in both the absence and presence of input correlations. The newly introduced sensitivity indices also quantify the importance of input independence and correlations, allowing one to determine whether or not the input correlations should be considered in practice. Quite recently the original Sobol's sensitivity indices are generalised to deal with correlated inputs, by using hierarchically orthogonal functional decomposition [START_REF] Chastaing | Generalized hoeffding-sobol decomposition for dependent variables-application to sensitivity analysis[END_REF][START_REF] Chastaing | Generalized sobol sensitivity indices for dependent variables: numerical methods[END_REF]. Compared to the generalised Sobol's indices, our indices are easier to understand and interpret.

Numerical examples and a practical application

In this section, analytic polynomial models, including one purely additive and three nonlinear ones, are taken as examples to illustrate the effectiveness and validation of our established analytic framework. A practical application of the method is also proposed to a deterministic HIV model. Ten involved parameters are then ranked according to their importance in establishing the uncertainty of the basic reproduction number R 0 .

Additive linear model

In the first example a purely additive model is investigated, with functional form given by

y = 2x 1 + x 2 + x 3 , (3.39)
where (x 1 , x 2 , x 3 ) ∼ N (µ, Σ) with mean vector µ = (0, 0, 0) T and covariance matrix which is constituted of the following fractional contributions of different dimensionality:

Σ =     1 ρ 12 2ρ 13
V 1 = 4, V 2 = 1, V 3 = 4, V 12 = 4ρ 12 , V 13 = 8ρ 13 , V 23 = 4ρ 23 , V 123 = 0. (3.42)
The vanishing nonlinear problem in the model function suggests the absence of coupling variance contributions but the presence of independent and correlated ones:

V U i = c 2 ii V i , V U j 12 = 0, V U k 13 = 0, V U l 23 = 0, V C i = (1 -c 2 ii )V i , V C j 12 = V 12 , V C k 13 = V 13 , V C l 23 = V 23 , V U i 123 = V C i 123 = 0, (3.43)
c 11 = (1 -ρ 2 23 ) -1/2 (1 -ρ 2 12 -ρ 2 13 -ρ 2 23 + 2ρ 12 ρ 13 ρ 23 ) 1/2 , (3.44)
c 22 = (1 -ρ 2 13 ) -1/2 (1 -ρ 2 12 -ρ 2 13 -ρ 2 23 + 2ρ 12 ρ 13 ρ 23 ) 1/2 , (3.45)
c 33 = (1 -ρ 2 12 ) -1/2 (1 -ρ 2 12 -ρ 2 13 -ρ 2 23 + 2ρ 12 ρ 13 ρ 23 ) 1/2 , (3.46)
in which {ρ 2 12 , ρ 2 13 , ρ 2 23 } ̸ = 1.

The underlying sensitivity measures are provided in table 3.1 under considering different correlations between input variables. Results indicate the vanishing sensitivity indices associated with the coupling effect between input independence and correlations. In the absence of correlated input variables (ρ = 0), the main sensitivity indices sum up to one.

By introducing input correlations, however, this summation could be smaller than one (with positive correlations) or larger than one (with negative correlations), contrary to the sum of the total sensitivity indices. A negative sensitivity index explains a negative partial variance contribution produced by the negative input correlation.

Nonlinear models

Trivariate model

In the second example a trivariate nonlinear model is considered, containing the linear, quadratic, and interaction terms: influential.

y = 2x 1 + x 2 2 + 4x 2 1 x 2 + x 1 x 3 , (3.47)

Fourvariate model

Another nonlinear model is designed based on four input variables as

y = x 1 x 3 + x 2 x 4 , (3.53)
where (x 1 , x 2 , x 3 , x 4 ) ∼ N (µ, Σ) with mean vector µ = (1, 2, 2, 1) T and covariance matrix which is constituted of In the calculation of V 1234 , the covariance cov(x 1 , x 2 , x 3 , x 4 ) of four variables is involved, whose derivation is presented in Appendix D. The form of model function (only involves the linear problem of each input) suggests the vanishing coupling effect in all partial variance contributions but the existent correlated and independent ones:

Σ =        1
V 1 = 4, V 2 = 1, V 3 = 1, V 4 = 4, V 12 = 4ρ
V U i = c 2 ii V i , V C i = (1 -c 2 ii )V i , V U j 13 = c 2 jj , V C j 13 = V 13 -V U j 13 , V U k 24 = c 2 kk , V C k 24 = V 24 -V U k 24 , (3.57)
where i ∈ {1, 2, 3, 4}, j ∈ {1, 3} and k ∈ {2, 4}. The partial variance contributions V 12 , V 14 , V 23 , V 34 and V 1234 are all contributed by input correlations. The coefficients c i are determined with Eq. (3.17). Table 3.

Ishigami function

The Ishigami function [START_REF] Ishigami | An importance quantification technique in uncertainty analysis for computer models[END_REF] has been extensively used as a benchmark for sensitivity analysis [START_REF] Morio | Global and local sensitivity analysis methods for a physical system[END_REF][START_REF] Lira | Beyond the gum: variance-based sensitivity analysis in metrology[END_REF]. Its functional form was defined as y = sin(x 1 ) + 7 sin 2 (x 2 ) + 0.1x 4 3 sin(x 1 ), (3.58) where all input variables are uniformly distributed in the interval [-π, π]. The presence of correlation between x 2 and each of the rest does not influence the total variance of In the first case, x 1 is formulated on the basis of x 3 as

x 1 = ρ 13 σ 1 σ 3 x 3 + √ 1 -ρ 2 13 r 1 , (3.59)
where the newly introduced random variable r 1 is the element of r that satisfies Eqs.

(3.14) and (3.15). Sensitivity measures show a strong positive variance contribution produced by the interaction effect between x 3 and the correlated part of x 1 , as well as a very strong negative variance contribution caused by the interaction term involving x 3 and both correlated and independent sections of x 1 .

For the second case, we generate x 3 on the basis of x 1 as

x 3 = ρ 13 σ 3 σ 1 x 1 + √ 1 -ρ 2 13 r 3 , (3.60)
where the random variable r 3 is the element of r that satisfies Eqs. (3.14) and (3.15).

Zero mean of x 1 leads to the nonexistence of sensitivity measures evaluating the main effect of x 3 . A dominated contribution to the variance of model response is produced by the interaction effect between x 1 and the coupling of independent section with the correlated one of x 3 .

Analysis results imply that the correlation between x 1 and x 3 , if exists, will play a crucial role in determining the model response. Sensitivity indices of correlated and that is considered in our method).

HIV model

The basic reproduction number, denoted as R 0 , is arguably the most important quantity in infectious disease epidemiology because it helps determine whether or not an infectious disease can spread through a population [START_REF] Fraser | Pandemic potential of a strain of influenza a (h1n1): early findings[END_REF][START_REF] Holme | The basic reproduction number as a predictor for epidemic outbreaks in temporal networks[END_REF]. R 0 is defined as the average number of new cases of an infection caused by one typical infected individual, in a population consisting of susceptibles only [START_REF] Anderson | Infectious diseases of humans: dynamics and control[END_REF][START_REF] Diekmann | On the definition and the computation of the basic reproduction ratio r 0 in models for infectious diseases in heterogeneous populations[END_REF][START_REF] Diekmann | The construction of nextgeneration matrices for compartmental epidemic models[END_REF]. The first application of this metric in epidemiology was introduced by George MacDonald in 1952, who designed agent-based models of the spread of malaria [START_REF] Macdonald | The epidemiology and control of malaria[END_REF]. Generally, the larger the value of R 0 , the harder it is to control the spreading of an epidemic. Typically, when R 0 < 1, the disease free equilibrium is locally asymptotically stable and the epidemic will die out in the long run, whereas if R 0 > 1, it is unstable and the epidemic will invade the population [START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF].

Consider a deterministic model of HIV-1 with vertical transmission (from an HIVinfected mother to her child) which was discussed in Ref. [START_REF] Usman | Sensitivity analysis of the parameters of a model of hiv-1 transmission[END_REF]. The basic reproduction number R 0 is represented by

R 0 = β 0 (1 -γ)θ 2 d + β 1 n 1 Q 0 (θ d -κ) + β 2 n 2 αQ 0 + (1 -γ)(κ + α)β 0 θ d θ d (θ d + κ)(θ d + α) . (3.61)
Description and baseline values of parameters included in the above expression are presented in table 3.5.

To identify the importance of individual parameters in establishing the uncertainty of R 0 , each parameter is artificially increased and decreased by 10% of its baseline value. The distribution of our analysis results is also presented in Fig. 3.5, along the direction of the correlation between β 1 and n 1 , and the direction of the correlation between β 2 and n 2 . The underlying results suggest that the effect contributed by the correlation

= 0.3, ρ β 2 n 2 = 0.5)
Table 3.7: The convergent analysis results of the basic reproduction number R0 for a deterministic HIV-1 model with vertical transmission by assuming independent and correlated input parameters. between β 1 and n 1 can be neglected, but the one contributed by the correlation between β 2 and n 2 is quite influential, in the estimation of the basic reproduction number R 0 .

x ρ = 0 ρ β 1 n 1 = 0.3, ρ β 2 n 2 = 0.5 sT i(= s U T i ) sT i s U T i s C T i Q0
Chapter 4

The establishment of sampling-based strategy

In diverse strategies for uncertainty and sensitivity analysis, the sampling-based method are considered as both effective and widely used in the absence of input correlations.

Analysis relying on the this method involves the generation and exploration of a mapping from input variables to the model response. The most crucial ingredient is the generation of sample points regarding sampling-based uncertainty and sensitivity analysis. In general, LHS and QMC sampling are often discussed and widely used for the analysis of computationally demanding models. The efficient stratification property of LHS allows for the extraction of a large amount of uncertainty and sensitivity information with a relatively small sample size [START_REF] Helton | Survey of sampling-based methods for uncertainty and sensitivity analysis[END_REF][START_REF] Helton | Latin hypercube sampling and the propagation of uncertainty in analyses of complex systems[END_REF]. Low-discrepancy sequences guarantee a perfect uniformity and uniqueness of QMC sampling in the input space, which avoids large computational cost in configuration average [START_REF] Kucherenko | Comparison of lagin hypercube and quasi monte carlo sampling techniques[END_REF][START_REF] Kucherenko | Estimation of global sensitivity indices for models with dependent variables[END_REF].

A comparison of our method with Sobol's one

A comparison of our method with Sobol's one is discussed in this chapter. Sobol's variance-based method is implemented here with the sampling-based strategy which also allows one to evaluate the sensitivity indices even if the functional form connecting input and output variables is absent. By using sample points, the representations of Sobol's sensitivity indices of different orders are presented in Section 1.3.7. In this part, the difference of our method from Sobol's one is presented via trivariate linear and nonlinear systems of specified functional relationships connecting y and x, in the absence of input correlations. Four numerical models of forms below are exemplified:

y 1 = 2x 1 -4x 2 + 6x 3 , (4.1)
y 2 = 2x 1 + 3x 2 2 + 4x 1 x 3 , (4.2)
y 3 = x 1 x 2 x 2 3 , (4.3)
y 4 = x 1 x 2 /x 3 , (4.4)
where input variables are assumed to be independent of each other and uniformly dis-

∞ ∑ i=1 1 i! ( ∂f ∂x 1 • ∂ i+1 f ∂x 1 ∂x i 3 )({µ}) • M 2 (x 1 )M i (x 3 ) / V (y), (4.5)
∞ ∑ i=1 1 i! ( ∂f ∂x 2 • ∂ i+1 f ∂x 2 ∂x i 3 )({µ}) • M 2 (x 2 )M i (x 3 ) / V (y), (4.6)
∞ ∑ i=1 1 i! ( ∂ 2 f ∂x 1 ∂x 2 • ∂ i+2 f ∂x 1 ∂x 2 ∂x i 3 )({µ}) • M 2 (x 1 )M 2 (x 2 )M i (x 3 ) / V (y). (4.7)
The above components are in sequence regarded as the effects of x 1 , x 2 and of the interaction effect between them in Sobol's definition since artificially moving the above items in our analytic formula for the third case yields {s 1 =0. E(y 1 )=2.00(00), E(y 2 )=3.00(00), E(y 3 )=0.08 [START_REF] Xu | Uncertainty and sensitivity analysis for models with correlated parameters[END_REF], E(y 4 )=0.86 [START_REF] Saltelli | Model-free importance indicators for dependent input[END_REF]; V (y 1 )=4.66 [START_REF] Williams | Safety evaluation and risk assessment of the herbicide roundup and its active ingredient, glyphosate, for humans[END_REF], V (y 2 )=2.57(66), V (y 3 )=0.01 [START_REF]Intergovernmental Panel on Climate Change. Climate Change 2014-Impacts, Adaptation and Vulnerability: Regional Aspects[END_REF], V (y 4 )=1.04 [START_REF] Marcot | Analysis of sensitivity and uncertainty in an individual-based model of a threatened wildlife species[END_REF].

as the square root of the whole information provided by x t alone which decides the sen-

sitivity index s t , that is √ V t .
Similarly, the contribution made by item [∂ 2 y/(∂x

1 ∂x 2 )]
in Eq. (4.7) can be approximated as √ V 12 . Nevertheless, the impact resulted from the right part in brackets of each equation among Eqs. (4.5), (4.6) and (4.7) should be separately approximated as half of the square root of V 13 , V 23 and V 123 . This because of the only existent even-order central moments of x 3 . x 3 is uniformly distributed, resulting non-existence of contributions from the relevant odd-order derivatives of f with respect to x 3 . Following these approximations, a modification is introduced to the Sobol's method:

s r 1 = s 1 - √ s 1 • s 13 2 , s r 2 = s 2 - √ s 2 • s 23 2 , s r 12 = s 12 - √ s 12 • s 123 2 , s r 13 = s 13 + √ s 1 • s 13 2 , s r 23 = s 23 + √ s 2 • s 23 2 , s r 123 = s 123 + √ s 12 • s 123 2 . (4.8)
The revised results for the last two cases are shown in y r 3 -and y r 4 -th rows in table 4.2. They are in well agreement with analytic values accompanied by ignorable errors. The above modification depends upon the functional form of the model under discussion that provides information about which of measures are overvalued and which are underestimated. Following the modification, it is easy to understand the difference between our method and Sobol's one. For models of absent nonlinear interaction effects of input parameters, our method (in the absence of input correlations) is equivalent to Sobol's one in evaluating the importance of individual input parameters. When the model under discussion involves nonlinear interaction terms, however, the difference between our method and Sobol's one is existed in quantifying some sensitivity indices. By introducing the modification to Sobol's method, a rough sampling-based method is also established, coincident with our analytic concept. This builds a foundation of extending the sampling-based method from independent case to the correlated one. Samplingbased method of input correlations is of great importance in practice, which allows one to decide where or not the input correlations should be considered in models of absent functional forms. Results also suggest that QMC sampling, compared to LHS, can provide more precise measures at lower computational cost, as stated in the first chapter. E(y 1 )=2.00(00), E(y 2 )=3.00(00), E(y 3 )=0.08 [START_REF] Xu | Uncertainty and sensitivity analysis for models with correlated parameters[END_REF], E(y 4 )=0.86 [START_REF] Most | Variance-based sensitivity analysis in the presence of correlated input variables[END_REF]; V (y 1 )=4.66(70), V (y 2 )=2.57 [START_REF] Orosz | Traffic jams: dynamics and control[END_REF], V (y 3 )=0.01(53), V (y 4 )=1.04(34).

Analysis of SIR and SIS models

The modeling of infectious diseases provides a tool to study the mechanisms by which diseases spread, to predict the future course of an outbreak and to evaluate strategies to control an epidemic [START_REF] Daley | Epidemic Modeling: An Introduction[END_REF]. Mathematical modeling of disease spreading was first carried out in 1766 by Daniel Bernoulli who designed a numerical model to defend the practice of inoculating against smallpox [START_REF] Bernoulli | Essai d'une nouvelle analyse de la mortalité causée par la petite vérole et des avantages de l'inoculation pour la prévenir[END_REF]. However, modern theoretical epidemiology did not begin until the research of Ronald Ross into the spread of malaria [START_REF] Ross | The prevention of malaria[END_REF]. The reliability investigation regarding deterministic compartmental models are performed for many years under the consideration of output's dependence upon one input factor at a time in fixing the others [START_REF] Wang | Epidemic spreading on complex networks with general degree and weight distributions[END_REF][START_REF] Wei | A unified framework of interplay between two spreading processes in multiplex networks[END_REF]. This widely used one-at-a-time method, however, fails to globally quantify the roles played by individual input factors and by their interaction effects in the estimation of model response. Following the deficiency, a systematical framework is carried out for the reliability analysis of both SIR and SIS models, by using sampling-based strategy and our proposed modification.

Considering a fixed population: N = S(t) + I(t) + S(t), McKendrick and Kermack designed the following nonlinear system of ordinary differential equations:

         ds(t) dt = -βs(t)i(t), di(t) dt = βs(t)i(t) -γi(t), dr (t) dt
= γi(t).

(4.9)

Apparently, we have s(t) + i(t) + r(t) = 1. β labels the probability at which a susceptible individual can be infected after communicating with infected individuals, and γ the probability at which an infected individual will be spontaneously recovered from the infection. In Eq. (4.9), dividing the first equation by the third one yields:

ln s(t) = - β γ r(t) + ln s(0) + r(0), (4.10)
in which s(0) and r(0) label the initial condition that specifies the infectious invading.

In general, r(0) is set to 0, yielding i(0) = 1 -s(0). At equilibrium state of disease spreading, the normalized number of both susceptible and infectious agents satisfy

         βsi = 0, βsi -γi = 0, γi = 0. (4.11)
which, together with Eq. (4.10), provides

   i = 0, ln s = -β γ (1 -s) + ln s(0). (4.12)
Apparently, SIR model does not lead to any analytic solution for the equilibrium state but just states that s associates with input parameters β, γ and s(0) which are assumed to be independent of each other and uniformly distributed in the real range [0, 1]:

s = f (s(0), β, γ); β, γ, s(0) ∈ [0, 1]. (4.13)
Differential equations for SIS model can be deduced from SIR model by removing the recovered populations into the susceptible compartment:

   ds(t) dt = -βs(t)i(t) + γi(t), di(t) dt = βs(t)i(t) -γi(t). (4.14)
Its exact solution at equilibrium state is Regarding SIS model yielding the exact solution of equilibrium state, the inaccuracy of the original sampling-based strategy in the implementation of our analytic method is introduced by the following summation set:

s =    1, γ ≥ β, γ β , γ < β.
∞ ∑ i=1 1 i! ( ∂s ∂γ • ∂ i+1 s ∂γ∂ i β )({µ}) • M 2 (γ)M i (β) / V (y), (4.16)
which leads to the overestimation of index s γ but underestimation of index s βγ . Based on this description, the relevant items are then corrected as

s r γ = s γ - √ s γ • s βγ 2 , s r βγ = s βγ + √ s γ • s βγ 2 . (4.17)
For SIR model with no exact solution of equilibrium state, the association of model response s with each input variable is shown in Fig. 4.2. It is noticed that s nonlinearly depends upon the ratio of γ to β but almost linearly correlates with s(0). Furthermore, the interaction effect of s(0) with each of β and γ can be neglected since the lines in double log plane of s-s(0) are almost parallel to each other for different ratios of γ to β except as s(0) → 1. In the region of s(0) → 1, interaction effects of s(0) with other factors are existed, explaining the non-vanishing high-order sensitivity indices involving s(0). Consequently, for SIR model, the items that give rise to the inaccuracy of the sampling-based strategy in the implementation of our analytic method include

∞ ∑ i,j,k=1 1 i! • j! • k! ( ∂ i s ∂γ i • ∂ j+k s ∂γ j ∂β k )({µ}) • M i+j (γ)M k (β) / V (y), (4.18)
∞ ∑ i,j,k=1

1 i! • j! • k! ( ∂ i s ∂β i • ∂ j+k s ∂γ j ∂β k )({µ}) • M j (γ)M i+k (β) / V (y), (4.19)
which yield the modification of sensitivity indices as Chapter 5

s r β = s β - √ s β • s βγ 2 , s r γ = s γ - √ s γ • s βγ 2 , s r βγ = s βγ + √ s β • s βγ 2 + √ s γ • s βγ 2 .
Opinion formation based on a gambling mechanism and its sensitivity analysis

Introduction

The dynamics formation of public opinion has raised a fair amount of attention from researchers working on employing diverse tools furnished by statistical physics to understand the appearance of some social phenomena [START_REF] Galam | Contrarian deterministic effects on opinion dynamics:the hung elections scenario[END_REF][START_REF] Barabasi | The origin of bursts and heavy tails in human dynamics[END_REF][START_REF] Ben-Naim | Opinion dynamics: rise and fall of political parties[END_REF]. As Castellano et al have stated, a common theme of social dynamics is the understanding of the transition from an initially disordered state to a configuration that displays order (at least partially) [START_REF] Castellano | Statistical physics of social dynamics[END_REF]. In opinion dynamics, ordered configuration corresponds to the emergence of an agreement which is necessary to be reached for a community in many situations appeared in daily life. In establishing a framework on opinion dynamics, the most important part includes specifying any possible opinion states of a community and defining the elementary processes or rules that determine opinion transitions of every individual between such states.

The first model of opinion dynamics known so far was proposed by Weidlich to study individual's decision behaviour by reformulating statistical models which are originally used for describing the dynamics of interacting spins [START_REF] Weidlich | The statistical description of polarization phenomena in society[END_REF]. Later on, Ising model made its first appearance in opinion dynamics [START_REF] Galam | Sociophysics: A new approach of sociological collective behaviour. i. mean-behaviour description of a strike[END_REF]. In these early models, as well as many other subsequently built ones, such as the widely considered majority-voter model [START_REF] Zhu | Opinion evolution on a ba scaling network[END_REF][START_REF] Pereira | Majority-vote model on random graphs[END_REF], and the Sznajd model [START_REF] Sznajd-Weron | Opinion evolution in closed community[END_REF], the public opinion is simplified as a discrete variable.

Discrete opinion often takes two possible values, representing a reasonable description in several instances: yes or no, buying or selling, cooperation or defection, Samsung or iPhone, etc. The discrete opinion dynamics is favoured over continuous one where agent's opinion can vary smoothly from one extreme of the range of possible choices to another, because of its simplicity in theoretical analysis [START_REF] Guo | Heterogeneity and subjectivity in binary-state opinion formation systems[END_REF][START_REF] Krapivsky | Dynamics of majority rule in two-state interacting spin systems[END_REF][START_REF] Sood | Voter model on heterogeneous graphs[END_REF]. However, for the investigation of the evolution of agent's attitudes about a given topic, for example, the political orientation of individual agents, continuous opinion dynamics seems more suitable. Furthermore, the latter is also definitely advantageous for gaining insight into the emergence of opinion clusters.

Opinion clusters could be one (consensus), two (polarization), or more (fragmentation).

How social clusters emerge in opinion space is one of the fundamental problems in social science. In small communities, such as a football team, or a research group, each agent knows and can communicate with every other one. For a large community or population, however, many agents do not know or even never heard of each other because of the geographic limitation. Hence, local interactions happen only within groups of friends. In computer simulations, the friendship in a community is defined in general by a connected network with the same size of the community. Moreover, the bounded confidence is also widely considered. It was proposed according to the realistic aspect of human communications happening among agents whose opinions are sufficiently close to each other [START_REF] Deffuant | Mixing beliefs among interacting agents[END_REF][START_REF] Hegselmann | Opinion dynamics and bounded confidence models, analysis, and simulation[END_REF][START_REF] Lorenz | Continuous opinion dynamics under bounded confidence: A survey[END_REF]. Mathematical models and statistical physics have offered powerful tools for deeply understanding how and what conditions collective behaviours happen in local individual interaction rules [START_REF] Degroot | Reaching a consensus[END_REF][START_REF] Chatterjee | Towards consensus: Some convergence theorems on repeated averaging[END_REF]. The roles of some possible social factors, such as leaders and environmental noises [START_REF] Zhao | Bounded confidence opinion dynamics with opinion leaders and environmental noises[END_REF], mass media [START_REF] Boudin | Opinion dynamics: Kinetic modelling with mass media, application to the scottish independence referendum[END_REF], agent's conviction and influence strength [START_REF] Vieira | Consequences of nonconformist behaviors in a continuous opinion model[END_REF], and informed agents [START_REF] Fan | Opinion evolution influenced by informed agents[END_REF], are also identified recently in the formation of opinion dynamics.

In current social opinion research, unidirectional influence from considered social factors to the underlying dynamics is the main focus. However, agent's opinion may also, in turn, react on the associated social environment. Based on the consideration, a virtual gambling mechanism is introduced to explore the characteristics of opinion dynamics behaviour in a social community. At each round of gambling, agents fight for a limited but conserved resource. The social pattern considered here is the resource distribution in the underlying system. Opinions of agents determine directly the resource that each agent can share, thereby deciding winning or losing of agents at our gambling. Agents then can accumulate scores from a reward and punishment system which is proposed based on the winning or losing state of agents. Agent's scores will, in turn, affect his opinion in subsequent interactions with his friends whose opinions must differ from his own not more than a pre-defined confidence threshold. In general, human beings are recognised as "higher animal" because of self-awareness and the freedom from nature's determinism that allows one to choose, whether for good or ill. This implies that agents participating in our gambling will instinctively learn from their friends who have earned more scores from the activity owing to the primary intend of winning money and/or material goods. A systematic framework is also performed for the uncertainty and sensitivity analysis of the model under discussion. The underlying analysis results help gain insight into how involved parameters can be ranked according to their importance in establishing the uncertainty of model response of our interest.

Modeling

In our model, we have N social agents. Their friendship is realized on a Barabási-Albert scale-free network (for detailed generation process, please check Ref. [START_REF] Albert | Statistical mechanics of complex networks[END_REF]) as many social networks have been reported to be scale-free [START_REF] Albert | Statistical mechanics of complex networks[END_REF][START_REF] Onnela | Structure and tie strengths in mobile communication networks[END_REF][START_REF] Choromański | Scale-free graph with preferential attachment and evolving internal vertex structure[END_REF]. Each agent i has an opinion s i in the real range [-1, 1]. Agents are segregated into two cliques according to the sign of their opinions: G + (with positive opinions) and G -(with negative opinions).

The absolute value of s i is regarded as the degree of attention or the cost that agent i pays for the issue of gambling.

A parameter γ is introduced by the gambling mechanism to control the resource discrepancy between two cliques: γ = R + /R -in keeping total resource R (=R + + R -) conserved. R + and R -label the resource allocated to G + and G -, respectively. In computer simulations, they are specified by γ and R defined previously:

R + = γR 1 + γ , R -= R 1 + γ , (5.1)
where the total resource R is set to 1 in our simulations, but the behavior of our model is almost independent of the value of R. For simplicity but without loss of generality, γ is restricted to the closed interval [0, 1] (the system is symmetric around γ = 1). The resource held by each clique is then shared by their members. Generally speaking, an agent who pays more attention or cost will share more resource according to the general rule of survival: the more you give, the more you receive. On this basis, the resource that an agent shares from his community is proportional to the absolute value of his opinion.

Rewards and punishments system is also introduced. An agent is regarded as a winner and to be rewarded with C scores if the resource he shares is no less than the global average value A r = R/N ; otherwise, the agent is considered as a loser and to be punished with losing C scores. C is a constant for each agent and set to 5 in our modelling, but its value does not influence the results of our interest. Actually, similar ideas have been considered in earlier modeling of social behaviors. In the naming game, as an example, the concept of score is introduced to indicate the agent's reputation which is variable in time and partially controls the behavior of the system [START_REF] Brigatti | Consequence of reputation in an open-ended naming game[END_REF]. For another example, in the modeling of discrete opinion dynamics based on a majority rule, a parameter indicating the convincing power of each agent, which may increase or decrease with time, analogous to the above mentioned score, is proposed to decide whether the majority opinion succeeds or not in a randomly selected group [START_REF] Crokidakis | The first shall be last: Selection-driven minority becomes majority[END_REF]. More considerations of the concept or similar concept of score are given by [START_REF] Nowak | Evolution of indirect reciprocity by image scoring[END_REF][START_REF] Crokidakis | Consequence of reputation in the sznajd consensus model[END_REF][START_REF] Crokidakis | The influence of local majority opinions on the dynamics of the sznajd model[END_REF].

In general, gamblers are economic with making decisions by their own and taking the strategy in their favour. Some people, however, may act as mindless while facing some events or subjects possibly because the corresponding events or subjects do not interest them or they do not care about personal gain. Following the consideration, a parameter p is introduced to control the fraction of mindless agents included in the system under consideration.

The initial state of the system is assumed to be fully disordered, that is, at the beginning of the dynamics, each individual has an opinion drawn from the uniform distribution in the range [-1, 1]. Mindless agents are selected at random from the system with probability p. At a given time step t, the following microscopic rules control the formation of opinion dynamics:

(1) Agent i is mindless. In this case, the agent takes the average opinion of those friends whose opinions differ from his own no larger than a certain confidence threshold ε to form his opinion for the next time, that is

s i (t + 1) = δ |F (i,t)| s i (t) + (1 -δ |F (i,t)| ) 1 |F (i, t)| ∑ j∈F (i,t) s j (t), (5.2)
where | • | for a finite set denotes the number of elements, δ x = 1(0) for x = 0(> 0) the Kronecker delta function. The set F (i, t), which denotes the effective friends of agent i, is defined by

F (i, t) = { 1 ≤ j ≤ N, j ̸ = i| a ij = 1 and |s i (t) -s j (t)| ≤ ε } , (5.3)
in which | • | denotes the absolute value of a real number, a ij the element in adjacency matrix with a ij = 1(0) if there is an edge connecting i and j (otherwise). In the present paper, the confidence threshold ε is taken as constant in time and across the whole population, as discussed in Ref. [START_REF] Weisbuch | Interacting agents and continuous opinions dynamics[END_REF] where the heterogeneity of threshold and the adaptive threshold were also introduced, respectively, yielding the appearance of new features of opinion dynamics.

(2) Agent i is economic. i will update his opinion by taking into account the scores of himself and of his effective friends. s i keeps unchanged if i gains the highest scores compared with all j ∈ F (i, t); otherwise, s i is updated by the opinion (or average opinion) of the agent(s) who own(s) the highest scores in the set F (i, t). The rule for the evolution of s i is given by

s i (t + 1) = δ c i (t)m i (t) s i (t) + (1 -δ c i (t)m i (t) ) 1 |M (i, t)| ∑ j∈M (i,t) s j (t), (5.4)
with

m i (t) = max ( c i (t), {c j (t), j ∈ F (i, t)} ) (5.5)
and

M (i, t) = { j ∈ F (i, t)| c j (t) = m i (t) } , (5.6)
where δ xy = 1(0) if x = y(otherwise), and max() is a function to return the maximum value of a set of elements. c i (t) labels the scores that i has gained until time t, which will be refreshed accordingly by

c i (t + 1) = { c i (t) + 5, r i (t) ≥ A r , c i (t) -5, r i (t) < A r .
(5.7) r i (t) is the resource that i shares at time t, given by

r i (t) =    s i (t) T Ω + (t) • R + , s i (t) > 0, s i (t) T Ω -(t) • R -, s i (t) ≤ 0, (5.8)
in which T Ω + (t) and T Ω -(t) are the sum of positive and negative opinions, respectively,

T Ω + (t) = N ∑ i=1,s i (t)>0 s i (t), T Ω -(t) = N ∑ i=1,s i (t)≤0
s i (t).

(5.9)

In our modelling, the unit time (t = 1) is defined as N , which corresponds to a single update of each agent. Network parameters are set as N = 1000 and k a = 4. k a indicates the average degree of each agent, that is the average number of friends.

Methods

Time evolution of the underlying opinion dynamics is firstly presented in Fig. 5.1. The absolute global average opinion and local average opinions in both cliques are represented as where N + (N -) indicates the number of agents in clique G + (G -). Numerical results demonstrate the presence of steady state on the opinion evolution process.

Ω(t) = 1 N N ∑ i=1 s i (t) , Ω + (t) = T Ω + (t) N + (t) , Ω -(t) = T Ω -(t) N -(t) , (5.10)
To characterise the coherence of the collective state of the population, we employ

Ω s = ⟨ 1 N N ∑ i=1 s i ⟩ , (5.11)
where ⟨• • • ⟩ denotes average over configurations, computed at steady states. Ω s is known to be an order parameter which plays the role of the magnetisation in magnetic systems.

It is sensitive to the unbalance between positive and negative opinions. A collective state with a significantly nonvanishing value of Ω s indicates a symmetry-broken distribution of opinions. A very small value of Ω s (≃ 0) means a symmetric distribution of opinions.

And an exact value Ω s = 0 corresponds to an absorbing state where all agents share the neutral opinion s = 0 [START_REF] Anteneodo | Symmetry breaking by heating in a continuous opinion model[END_REF].

In practical simulations, first 100 time steps are used to make the system reach its steady state. Time average is then taken over the following 100 steps. 1000 independent simulations are proposed to take the average over configurations.

Results and discussion

The global steady state is driven from the local steady states of both cliques. This implies that the order parameter can be represented as

Ω s = n + * Ω + s + (1 -n + ) * Ω - s , (5.12)
where n + indicates the normalized number of agents in clique G + ; Ω + s and Ω - s denote the average opinions in clique G + and G -at steady state, respectively. Firstly, the role of resource allocation parameter γ in the formation of opinion dynamics is focused by fixing both confidence threshold and fraction of mindless agents. For simplicity, we set ε = 2 and p = 0. The primary intent of winning money or material goods at gambling will drive agents to learn opinions (or strategies) of their friends who have gained more scores. At the beginning of the dynamics (t = 0), rewards and punishments system requires the opinions of winners (have the ability to affect others) to follow which guarantees the presence of winners in G + at the initial time. We divide the opinion space into three adjacent regions, as shown in Fig. 5.2. An agent i with initial opinion s i (0) in the region R2 will adapt his opinion to R1 or R3, or continuously to R2 by learning from his friends. However, to keep s i (t > 0) staying in R2, i must have friends with initial opinions in neither or both of the regions R1 and R3. This event is neglected in the following analysis as it occurs with a quite small probability decided by the power function of the fraction of agents in different regions.

|s i (0)| ≥ { 1+γ 4γ , s i (0) > 0,
Fully disordered initial state also drives uniformly distributed opinions in both regions R1 and R3 at very beginning of gambling. With agents migrating in opinions from R2

to R1 or R3, the resource shared by each agent will decrease accordingly. This gives rise to a higher requirement for the opinions of agents who can continuously win at our gambling. The final range of opinions in both cliques can be approximated by

s i N + • 1 2 (1 + 1+γ 4γ ) • R + ≥ R N , (5.15)
contribution of other agents is negligible. This provides

Ω s = 1 1 + γ .
(5.22)

When we limit γ to [10/N, 1/3), the fraction of winers is close to one. This implies two main clusterings of agents in opinions -1 and 1 because the condition of winning requires the final minimal absolute opinions in both cliques to satisfy

s N + • 1 2 (1 + s) • R + ≥ R N , for s > 0, (5.23
)

and |s| N -• 1 2 (1 + |s|) • R -≥ R N , for s < 0, (5.24)
respectively. Clustering of few agents is ignored in opinion close to 0 here, which is generated by agents who have friends in both above mentioned clusters. The order parameter Ω s is then specified as

Ω s = | γ 1 + γ - 1 1 + γ | = 2 1 + γ - 1.
(5.25)

Now we analytically obtain the distribution of order parameter Ω s in the whole region of resource allocation parameter γ by keeping p = 0 and ε = 2:

Ω s =              1, 0 ≤ γ < 1 N , 1 1+γ , 1 N ≤ γ < 10 N , 2 1+γ -1, 10 N ≤ γ < 1 3 , 3 2(1+γ) -3 4 , γ ≥ 1 3 .
(5.26)

In the present paper, the size of system considered is N = 1000. This implies the above equation can be simplified as

Ω s = { 2 1+γ -1, 0 ≤ γ < 1 3 , 3 2(1+γ) -3 4 , 1 3 < γ ≤ 1.
(5.27)

A remarkable agreement between analytic values and numerical results is presented in p = γ = 0

The effect of confidence threshold ε is also assessed on the formation of opinion dynamics by assuming p = γ = 0. γ = 0 stipulates the total resource is concentrated in one clique (G -). Following Eq. (5.13), initial winners hold opinions between -1 and -0.25. By introducing confidence threshold, an agent i with initial opinion larger than -0.25 can be affected at a probability P i by his friends with opinions between -(ε -s i ) and -0.25.

P i is straightforwardly specified by the mean filed assumption as

P i = k a 2 (ε -s i (0) - 1 4
), (5.28) where k a (=4) denotes the average degree. Setting P i = 1 yields

s i (0) = ε - 1 4 - 2 k a = ε - 3 4 .
(5.29)

We divide the opinion space into four adjacent regions, as shown in Fig. 5.5. Agents with initial opinions in region R2 will definitely adapt their opinions to R1 by local communications. Agents who take initial opinions in R3 will be affected by their effective friends at a certain probability. For the fourth region, opinion distribution will maintain is approximated as

n 3→1 = ∫ ε-1 4 ε-3 4 1 2 • k a 2 (ε -x - 1 4
)dx = 1 8 .

(5.30)

The fraction of winners (all agents in G -) in the steady state is then determined by

n w = 1 8 + 1 2 (ε - 3 4 + 1) = ε 2 + 1 4 , (5.31)
which says the final opinions of agents in G -must satisfy the inequality

|s| N ( ε 2 + 1 4 ) • 1 2 (|s| + 1) • R -≥ R N , (5.32)
where opinions in G -are assumed to be uniformly distributed between -1 and s. Sub-

stituting R -= R drives s ≤ 2ε + 1 2ε -7 , (5.33)
which also gives the final average opinion in G -:

Ω - s = 2ε + 1 2ε -7 .
(5.34)

The final average opinion in R4 is readily calculated:

Ω 4 s = 1 2 [1 -(ε - 1 4 )] • 1 2 [1 + (ε + 1 4 )] = - ε 2 4 + ε 8 + 15 64 .
(5.35)

Regarding region R3, some agents will be affected by his friends with initial opinions in R1. Furthermore, agents with opinions more closer to R2 are more easier to be affected. For fully connected networks, we have

Ω s =        16 13-4ε [ 1 + 1 4 (ε -1 4 ) 2 ] -5 4 , 0 ≤ ε ≤ 1 4 ; 16 13-4ε + 1 4 (ε -1 4 ) 2 -5 4 , 1 4 < ε ≤ 5 4 ; 1, 5 4 < ε ≤ 2.
(5.41)

A remarkable agreement between analytic values and numerical results is presented in Fig. 5.6, attached with a reasonable discrepancy when the friendship network is sparse (Fig. 5.6(a)). The small discrepancy is caused by mean field assumption, which gradually decreases with network parameter k a (check plots from left to right) and will vanish when the underlying network is fully connected, i.e., when the mean field limit reaches.

The distribution of Ω s in general situations

Consider different values of both p and γ. We would like to investigate the effect of parameter ε on the formation of opinion dynamics in detail. The underlying numerical results are exhibited in Fig. 5.7: the first panel relates to the distribution lines for different p by keeping γ = 0, and the second one involves the data for different γ by setting p = 0. For any fixed p and γ, Ω s is positively correlated with ε. This is because higher confidence threshold guarantees the happening of more local communications between an arbitrary agent and his friends in the system under consideration, which also benefits the formation of more ordered opinion dynamics. Equation (5.13) suggests a vanishing cross-opinion influence between two cliques when ε ≤ (1 + γ)/4. Consequently, the fraction of agents whose opinions will be affected by others can be approximated as

n f = ∫ ε-1+γ 4 -1+γ 4 1 2 • k a 2 (ε -x - 1 + γ 4 )dx + ∫ 1+γ 4γ 1+γ 4 -ε 1 2 • k a 2 (x + ε - 1 + γ 4γ )dx = k a 4 ε 2 .
(5.42)

Apparently, for sparse networks (k a is small), the evolution of opinion dynamics is neg- between before and after introducing mindless agents by considering both low and high confidence thresholds. p = 0 p = 0 says all agents involved in our gambling are economic. Numerical results suggest the presence of platform in the distribution of Ω s for fixed γ (Fig. 5.7(b)). This implies, when ε is sufficiently high, its further increase will not influence the established dynamics.

Denote the critical confidence threshold as ε c from which a platform is generated.

Reviewing the special case with γ = 0, agents will try to migrate in opinion from G + to G -for gaining scores. Eq.(5.29) states that the initial opinion of an agent i from G + who will be definitely affected by others must follow

s i (0) ≤ ε - 1 4 - 2 k a . ( 5

.43)

To ensure the migration of all agents in opinion from G + to G -, we have

1 ≤ ε - 1 4 - 2 k a , (5.44)
which determines

ε c = { 5 4 + 2 ka , k a ≥ 8 3 , 2, k a < 8 3 .
( .46) where k a = 4 is inserted. To ensure agent i can be definitely affected by others that hold opinions in [-1, -(1 + γ)/4], its initial opinion must satisfy

s i (0) = ε - 3 + γ 4 . ( 5
s i (0) ≤ ε - 3 + γ 4 .
(5.47)

By ignoring the migration of agents with initial opinions larger than ε -(3 + γ)/4 (rare event), the critical confidence threshold then can be approximated as

1 2 (ε c - 3 + γ 4 ) = 1 2 - γ 1 + γ , (5.48)
where n + = γ/(1 + γ) (final fraction of agents in G + ) was used. Now we get

ε c = 3 + γ 4 + 1 1 + γ . (5.49)
For comparable resource allocation between two cliques: γ ∈ (1/3, 1], winners are produced in both cliques at the beginning of dynamics: s i (0) ≤ -(1 + γ)/4 or s i (0) ≥

(1 + γ)/(4γ). The average probability at which an agent i with initial opinion between -(1 + γ)/4 and (1 + γ)/(4γ) will be affected is represented by

P i = 1 2 [ k a 2 (ε -s i (0) - 1 + γ 4 ) + k a 2 (ε + s i (0) - 1 + γ 4γ ) ] . (5.50)
Inserting k a = 4 and setting P i = 1 provide

ε c = 1 2 + (1 + γ) 2 8γ .
(5.51)

Opinion Clustering

A dispersion index proposed by Derrida and Flyvberg [START_REF] Derrida | Multivalley structure in kauffman's model: Analogy with spin glasses[END_REF] is introduced to identify agents clustering in opinions: It is noticed that, in the absence of economic agents, Y is a fast increasing function of parameter ε. The emergence of large cluster compared to which other clusters, if exist, are negligible (R(Y, m L ) ≈ 0) happens from ε = 0.4. Furthermore, the consensus state is definitely generated since ε = 1. When economic agents are present, however, the emergence of consensus state requires higher confidence threshold ε and p = γ = 0. Introducing of mindless agents or increase of the difference in resource between two cliques will destroy the clustering of agents in one opinion so as to generate small clusters with comparable size at high ε. Nevertheless, when ε is not so large (< 1.8), dispersing resource to two cliques is better for opinion clustering, compared to concentrating resource in one clique where comparable clusters in size are produced in adjacent opinions.

Y = M ∑ j=1 m 2 j , (5.52
Colour fill diagrams of dispersion index Y (Fig. 5.10) present an intuitive interpretation for the roles of both p and γ in the clustering of agents in opinions. For fixed ε, opinion clustering happens only in the field with small p and γ (p < 0.5, γ < 0.4) and when p tends to 1. However, this clustering will fade away with the decrease of ε. Along with p-axis, the smallest value of Y is obtained at p ≈ 0.7 which corresponds to the most chaotic steady state of opinion dynamics.

The fraction of winners

Regarding gambling game, one of most attractive observations is the fraction of winners: n w . Numerical results of n w are presented in γ-ε space by fixing parameter p, as shown in Fig. 5.11(a). In the upper half region of ε, increase of mindless agents will decrease the generation of winners. However, when ε is fixed between 0.5 and 1, a fragile positive correlation between n w and p exists at low γ. Furthermore, for small p, n w acts as a positive function of ε for each specified γ except the pattern formed at the position with ε ≈ 1.5 and γ ≈ 0.1. The pattern, corresponding to the presence of a small number of winners, will gradually disappear with the increase of p, accompanied by the appearance of a non-monotonic relationship between n w and ε. Robust association of n w with parameter γ happens only when p is large and ε > 0.5. A weak positive correlation between n w and γ appears in the region of low confidence threshold (ε < 0.5) (opinions of almost all agents keep their initial states), which is introduced by fully disordered initial dynamics that gives

n w (0) = { 3-γ 8 , γ ≤ 1 3 , 1 -(1+γ) 2 8γ , γ > 1 3 .
(5.54)

The association of the fraction of economic agents in the group of winners, denoted by N e w , with that in the whole system (= 1 -p) is also investigated for identifying who, an economic one or a mindless one, is more likely to win at our gambling game (Fig.

5.11(b))

. Results display that ε plays a fragile role in deciding the probability of winning for economic or mindless agents. The increase of γ, however, enables a general economic agent to win with higher probability. Furthermore, values of n e w are located only in the upside of diagonal (n e w = 1 -p, corresponding to the equal opportunity of winning for economic and mindless agents). This suggests, in general case, it is more likely for economic agents than mindless ones to win at gambling.

Uncertainty and Sensitivity analysis

Finally, a systematic framework on sensitivity analysis of the model under discussion is performed by the use of sampling-based strategies including LHS and QMC sampling. Particularly, QMC sampling is implemented by using the high-dimensional Sobol low-discrepancy sequence generator. Sobol low-discrepancy sequence holds advanced uniformity properties in the input space, and makes QMC sampling deterministic for previously defined input space, as stated in the first chapter. This avoids highly computational cost in performing configuration average. Regarding nonlinear models, sensitivity analysis allows one to globally understand the importance of individual input factors and of their interaction effects in establishing the uncertainty of model response. The model response of our interest here includes the order parameter Ω s of opinion dynamics, the dispersion index Y of opinion clustering, and the fraction n w of winners at gambling. Analysis results are listed in table 5.1. LHS-based the resource allocation parameter γ, confidence threshold ε, and the interaction effect between them. In establishing the uncertainty of dispersion index Y , parameter p, the fraction of mindless agents, acts as the most important factor. Furthermore, the interaction effects among three involved parameters are non-ignorable in determining the uncertainty in Y . Parameter p also plays the most important role in the generation of winners at our gambling, compared to which, the influence contributed by each of the remaining parameters alone is ignorable.

Chapter 6

Conclusions and Prospectives

Conclusions

The present research focuses on the establishment of a systematic theoretical framework on the uncertainty and sensitivity analysis of complex systems, as well as its applications to some numerical models. Specifically, the remarkable results and potential applications are concluded as below.

Firstly, an analytic formula is deduced by using Taylor series to exactly evaluate the uncertainty in output variable with given uncertainties in input factors in the absence of input correlations. With its applications to univariate power-law and exponential functions, it is analytically demonstrated that the widely used approximation with only the first order contribution of input uncertainty considered can satisfactorily explain the output variance only when the input uncertainty is negligible or the input-output relationship is close to linear. This implies the importance of high-order contributions of input uncertainty in the analysis of highly nonlinear models. Applications of the analytic method are proposed to the power grid system and EOQ model where input factors are assumed to be statistically independent of each other. Input factors of both models are analytically ranked according to their importance in establishing the uncertainty of model response.

The analytic method is then extended to the situation of input correlations. An arbitrary variable can be divided into independent and correlated parts. Coefficients that specify the correlated and independent parts are then derived by using linear correlation model. With given correlation coefficients between any two input parameters, it is straightforward to quantify the sensitivity of model response with respect to the input independence and correlations, as well as to their coupling effects. This allows one to decide whether or not the input correlations should be considered in practice. Four numerical examples have confirmed the effectiveness and applicability of the extended analytic method. Specifically, the method is applied to a deterministic HIV model.

Analysis results provide the importance of ten involved factors in determining the basic reproduction number R 0 . This may help effectively limit the spread of HIV by controlling the first three most important parameters: the rate of development to AIDS (κ), the contact rate of susceptible with symptomatic infective (β 2 ), and the number of sexual partners of susceptible with symptomatic infective (n 2 ), and simplify an HIV model by neglecting the effects of two negligible parameters: the birth rate of infective (β 0 ) and the fraction of susceptible newborn from infective class (γ). Moreover, the correlation between β 2 and n 2 , if exists, also provides a non-negligible effect on the basic reproduction number.

The difference of our method from Sobol's one is also discussed in the absence of input correlations. Sobol's method is implemented with the sampling-based strategy. Results

suggest that our method is equivalent to Sobol's one when the nonlinear interaction effects between input factors are absent in the model under discussion. However, when the model involves nonlinear interaction effects between input factors, Sobol's method may overvalue the roles acted by individual factors but underestimate those played by their interaction effects, compared to our analytic results. Following the phenomenon, a modification is introduced to Sobol's method for establishing a sampling-based strategy that well coincides with our analytic method. The modification depends on the functional form connecting input and output variables that tells which of sensitivity indices are overvalued and which are underestimated. Furthermore, both SIR and SIS models regarding the epidemic spreading dynamics are analysed by the modified sampling-based strategy. Results suggest that the equilibrium state of SIR model more robustly depends upon the initial proportion of susceptible individuals than both infectious probability β and recovered rate γ; their interaction effects totally contribute 30% of the information for the decision of equilibrium state. As for SIS model, parameter γ plays more crucial role than β, by explaining 49% variance in observation s. Their interaction effect also provides non-ignorable (31%) contribution to the determination of epidemic spreading.

Finally, we design a model based on a virtual gambling mechanism to investigate the dynamics formation of public opinion in the presence of limited but conserved resource and of confidence threshold. Theoretical analysis based on the traditional one-at-a-time method provides a deep understanding of the roles of both resource allocation parameter and confidence threshold in the formation of public opinion dynamics. For a sparse network, the evolution of opinion dynamics is negligible in the period of low confidence threshold when mindless agents are absent. Furthermore, the consensus state which also means a win-win state is generated only when the following three conditions are satisfied simultaneously: mindless agents are absent, the resource is concentrated in one clique, and confidence threshold tends to a critical value which is negatively dependent upon the average number of friends of individual agents. Numerical results also imply that, for specified fraction of mindless agents, the dependence of opinion dynamics upon confidence threshold is weakened with the decrease in resource difference between two cliques. The clustering of agents in opinions is also considered. The introducing of gambling mechanism leads to a higher requirement for confidence threshold (> 1.2) to drive apparent clustering of agents, which also demands the percentage of mindless agents to be less than half and resource allocation between two cliques to be very biased. For the low fraction of mindless agents (<0.5), resource partitioning is more advantageous to opinion clustering when confidence threshold departs from 2, compared to resource concentration. In addition, for fixed γ and ε, the most chaotic steady state of opinion dynamics is generated at p ≈ 0.7. It is also stated that, in general cases, it is more likely for economic agents to win at gambling, instead of mindless ones. By the use of sampling-based methods, a systematic framework for the uncertainty and sensitivity analysis of the model under discussion is provided. Results state that the resource allocation parameter, confidence threshold, and their interaction effect all play robust roles in establishing the uncertainty of the order parameter of opinion dynamics. Regarding the dispersion index of opinion clustering and the fraction of winners, their uncertainties are mainly contributed by the fraction of mindless agents, and also by the interaction effects among three involved parameters. This research may contribute to the deep understanding of the generation of some social dynamics, and of the roles played by any parameters that drive the underlying dynamics.

Future work plan

Following the above research work for uncertainty and sensitivity analysis of complex systems, several aspects may deserve further considerations, whereby the roles of input parameters can be more precisely quantified in both the modeling of complex systems and the analysis of practical problems.

Sampling-based strategy design in the presence of input correlations

A theoretical framework has been established for the uncertainty and sensitivity analysis of general complex systems with given functional relationships connecting input and output variables. However, in the modeling of most complex systems (e.g., the agentbased system), it is often hard to specify the exact functional form of models. In the absence of functional forms connecting input and output variables, the uncertainty and sensitivity analysis should be performed only through sampling-based strategies. At As a consequence, x 1 , if only correlated with x 2 , can be formulated as

x 1 = σ 1 σ 2 ρ 12 x 2 + √ 1 -ρ 2 12 r 1 . (B.6)

B.2 Three correlated variables

If correlations just exist among three variables, say x 1 , x 2 and x 3 , x 1 , for example, then can be similarly formulated by the use of linear correlation model as

x 1 = a 12 x 2 + a 13 x 3 + c 11 r 1 . (B.7)
Recalling the correlation coefficient between x 1 and x 2 , we get

ρ 12 = 1 σ 1 σ 2 E [(x 1 -µ 1 )(x 2 -µ 2 )] = 1 σ 1 σ 2 E [(a 12 (x 2 -µ 2 ) + a 13 (x 3 -µ 3 ) + c 11 (r 1 -µ(r 1 )))(x 2 -µ 2 )] = σ 2 σ 1 a 12 + σ 3 σ 1 a 13 ρ 23 . (B.8)
Analogously, the correlation coefficient between x 1 and x 3 can be expressed as

ρ 13 = 1 σ 1 σ 3 E [(x 1 -µ 1 )(x 3 -µ 3 )] = 1 σ 1 σ 3 E [(a 12 (x 2 -µ 2 ) + a 13 (x 3 -µ 3 ) + c 11 (r 1 -µ(r 1 )))(x 3 -µ 3 )] = σ 2 σ 1 a 12 ρ 23 + σ 3 σ 1 a 13 , (
σ 2 1 = E [ (x 1 -µ 1 ) 2 ] = E [ (a 12 (x 2 -µ 2 ) + a 13 (x 3 -µ 3 ) + c 11 (r 1 -µ(r 1 ))) 2 ] = a 2 12 σ 2 2 + a 2 13 σ 2 3 + 2a 12 a 13 σ 2 σ 3 ρ 23 + c 2 11 V (r 1 ). (B.11)
Inserting the expressions of a 12 and a 13 , and σ 2 1 = V (r 1 ), into the above equation, we obtain

c 11 = √ 1 -ρ 2 12 -ρ 2 13 -ρ 2 23 + 2ρ 12 ρ 13 ρ 23 1 -ρ 2 23 .
(B.12)

Accordingly, an arbitrary variable x 1 that is correlated with x 2 and x 3 at the same time can be generated with given pairwise correlation coefficients through equation

x 1 = ρ 12 -ρ 13 ρ 23 1 -ρ 2 23 σ 1 σ 2 x 2 + ρ 13 -ρ 12 ρ 23 1 -ρ 2 23 σ 1 σ 3 x 3 + √ 1 -ρ 2 12 -ρ 2 13 -ρ 2 23 + 2ρ 12 ρ 13 ρ 23 1 -ρ 2

B.3 Four correlated variables

Consider a more complex situation where correlations exist among four variables, say x 1 , x 2 , x 3 and x 4 . As before, employing the linear correlation model, we can represent

x 1 as x 1 = a 12 x 2 + a 13 x 3 + a 14 x 4 + c 11 r 1 , (B.14)
Reviewing the definition of correlation coefficient, we have where µ i = 0 for i = 1, 2, 3 were used. The above expression suggests a vanishing coupling variance contribution but the existent independent and correlated ones: where µ i = 0 for i = 1, 2, 3 were used. In the above equation, the first term is produced by the correlations of x 2 with the remaining inputs, the second one by its independence, and the third one by the coupling effect associated with correlations and independence, indicating

ρ 12 = 1 σ 1 σ 2 E [(x 1 -µ 1 )(x 2 -µ 2 )] , ρ 13 = 1 σ 1 σ 3 E [(x 1 -µ 1 )(x 3 -µ 3 )] , ρ 14 = 1 σ 1 σ 4 E [(x 1 -µ 1 )(x 4 -µ 4 )] . (
V U 1 = 4c
V U 2 = 2c 4 22 , V C 2 = 2(1 -c 2 22 ) 2 , V UC 2 = 4c

C.2 Second-order contributions

The second-order partial variance contributions produced by the combinations between each pair of inputs are expressed as ] , (C.12)

V 12 = 1 2! • 2! ( ∂ 3 y
V 13 = ( ∂ 2 y ∂x 1 ∂x 3 ) 2 ({µ}) • [ M 2,2 (x 2 1 , x 2 3 ) -M 2 1,1 (x 1 , x 3 ) ] = M 2,2 (x 2 1 , x 2
3 ) -ρ 2 13 , (C.13)

V 23 = 0. (C.14)
Considering the reforming of x 1 on the bases of x 2 and x 3 , the partial variance contribution V 12 is determined as Regarding the formulation of input x 2 , V 12 can also be calculated as Recalling the formulation of x 1 on the bases of x 2 and x 3 , the partial variance contribution V 13 can be analogously obtained as which only contains the correlated variance contribution (first three components) produced by the correlations of x 1 with the rest variables, and independent one (last component) produced by the independence of x 1 :

V 12 =16E [ x 2
V 12 = 16 [ x 4
V 13 = E [ x 2
V U 1 13 = c

C.3 Third-order contributions

The third-order partial variance contribution associated with the second test case is expressed as where µ i = 0 for i = 1, 2, 3 were used. In the above equation, the first two items are contributed by the correlations of x 2 with both x 1 and x 3 , and the last two items, summing to zero, are provided by the coupling effect between the independence of where (x 1 , x 2 , x 3 , x 4 ) ∼ N (µ, Σ) with mean vector µ = (1, 2, 2, 1) T and covariance matrix displayed in (3.54). Select x 1 (equivalent to select each of the others) to be reformed on the bases of the remaining factors: ) = { 2 i-1 (-1) (i+2)/2 , i is even, 0, i is odd. in the form of Ishigami function. x 2 is assumed to be independent. But x 1 and x 3 may correlate with each other. To interpret the independent, correlated and also coupling variance contributions included in V 1 , x 1 should be formulated on the basis of x 3 : ] , (E.9)

V 123 = 2 2! ( ∂ 2 y ∂x 2 2 • ∂ 2 y ∂x 1 ∂x 3 )({µ}) [ M 1,2,1 (x 1 , x 2 2 , x 3 ) -M 2 (x 2 )M 1,1 (x 1 , x 3 ) ] = M 1,
x 1 =
x 1 = ρ 13 σ 1 σ 3 x 3 + √ 1 -ρ 2 13
where µ 3 = µ(r 1 ) = 0 were used. By employing the binomial theory, the above expression is updated as .

V 1 = ∞ ∑ i,
(E.17)

V 13 is now obtained as The Independent and correlated partial variance contributions included in V 13 are specified by setting k = 0 and k = i + j + 1, respectively:

V U 1 13 = 0.1 ∞ ∑ i,j=0
(-1) i+j π 2(i+j+1)+4 (2i + 1)!(2j + 1)! (1 -ρ 2 13 ) i+j+1 2(i + j) + 3

[ 0. . (E.20)

The coupling variance contribution is then spontaneously determined for x 1 by

V UC 1 13 = V 13 -V U 1 13 -V C 1 13 . (E.21)
Optionally, x 3 can be also formulated in terms of x 1 as: . (E.25)

x 3 = ρ 13 σ 3 σ 1 x 1 + √ 1 -ρ 2 13
The independent and correlated effects contained in the above partial variance contribution can be calculated by setting k = l = 0 and k = 4, l = 2, respectively:

V U 3 13 = 0.1 ∞ ∑ i,j=0
(-1) i+j π 2(i+j+1)+4 (2i + 1)!(2j + 1)! [ 0. The coupling variance contribution is then determined for x 3 by

V UC 3 13 = V 13 -V U 3 13 -V C 3 13 .
(E.28) The primary intend of winning scores at gambling will impel agents with opinions in R1 to win. Denote the maximal steady opinion held by winners as s. The steady value of average opinion in G -then can be determined accordingly by

Ω - s = Ω 1 s + Ω 2 s + Ω 3 s = n 1 • 1 2 (s -1) - 1 6 ε 3 + 3 8 ε 2 - 1 8 ε - 1 4 (ε - 1 4 ) 2 .
(F.6)

The definition of winning at our gambling requires The order parameter is now directly stated as Ω s = 1 13 -4ε 2 ( 8 3 ε 3 + 3ε 2 ). (F.9)

When we limit the confidence threshold to the range (1/4, 1/2], the opinion space then can be partitioned as Fig. The stable average opinions in R2 and R3 are given by

Ω 2 s = ∫ 0 -1 4 1 2 [1 - k a 2 (ε -x - 1 4 )] • xdx = 1 32 (ε - 7 12
), (F.12)

and

Ω 3 s = ∫ ε-1 4 0 1 2 [1 - k a 2 (ε -x - 1 4 )] • xdx = - 1 6 (ε - 1 4 ) 2 (ε - 7 4
).

(F. [START_REF] Lee | Response surface augmented moment method for efficient reliability analysis[END_REF] respectively. In the region R4, opinions are fixed at fully disordered initial positions throughout the evolution of dynamics. This offers from R3 to R1 is approximated as

Ω 4 s = 1 4 [1 -(ε - 1
n 3→1 = ∫ 1 ε-3 4 1 2 • k a 2 (ε -x - 1 4 )dx = - 1 2 (ε - 5 4
) 2 + 1 8 .

(F.25)

Consequently, the fraction of agents in G -will fix at

n -= 1 2 (ε - 3 4 + 1) + n 3→1 = 1 - 1 2 (ε - 7 4 ) 2 .
(F.26)

The primary intend of winning scores at gambling will guide all agents in G -to win.

Analogously, the final average opinion in G -can be obtained as below by recalling Eq.

(F.7):

Ω - s = 1 - 4 (ε -7 4 ) 2 + 2 . (F.27)
For the clique G + , we have

Ω + s = ∫ 1 ε-3 4 1 2 [1 - k a 2 (ε -x - 1 4 )] • xdx = ε 6 (ε - 9 8
) 2 -55 128 (ε -1) + 5 24 .

(F.28)

The corresponding order parameter is now obtained as having initial opinions in R2 will migrate to R1 in opinion. Analogously, the agents with opinions in R1 then will learn from their friends to win in the opinion evolution process.

Ω s = 4 (ε -7 4 ) 2 + 2 - ε 6 (ε - 9

Abstract

By means of taylor series expansion, a general analytic formula is derived to characterise the uncertainty propagation from input variables to the model response, in assuming input independence. By using power-law and exponential functions, it is shown that the widely used approximation considering only the first order contribution of input uncertainty is sufficiently good only when the input uncertainty is negligible or the underlying model is almost linear. This method is then applied to a power grid system and the eoq model. The method is also extended to correlated case. With the extended method, it is straightforward to identify the importance of input correlations in the model response. This allows one to determine whether or not the input correlations should be considered in practical applications. Numerical examples suggest the effectiveness and validation of our method for general models, as well as specific ones such as the deterministic hiv model. The method is then compared to Sobol's one which is implemented with sampling based strategy. Results show that, compared to our method, it may overvalue the roles of individual input factors but underestimate those of their interaction effects when there are nonlinear coupling terms of input factors. A modification is then introduced, helping understand the difference between our method and Sobol's one. Finally, a numerical model is designed based on a virtual gambling mechanism, regarding the formation of opinion dynamics. Theoretical analysis is proposed by the use of one-at-a-time method. Sampling-based method provides a global analysis of output uncertainty and sensitivity.

Key Words

Uncertainty analysis, Sensitivity analysis, Variance decomposition, Sampling, Correlation, Sensitivity measure, Opinion dynamics L'Université Bretagne Loire

Investigation on uncertainty and sensitivity analysis of complex systems

Yueying ZHU

Figure 1 . 1 :

 11 Figure 1.1: Modeling after Rosen (1991)

Figure 1 . 2 :

 12 Figure 1.2: Uncertainty characterisation of independent input variables. (a): an example of given distribution laws of input variables: input 1 is uniformly distributed in the real range [0, 1]; input 2 follows the standard normal distribution (µ = 0, σ = 1). (b): an example of deterministic input variables: input 1 is specified at 1; input 2 is fixed at 0.5. Their uncertainties are represented by changing 50% around their normative values.

2 )

 2 It is straightforward to state the uncertainty in output y according to its values presented in Eq. (1.2).Presentation formats of the uncertainty in model output include mathematical expectation, standard deviation, the percentiles of its distribution, confidence bounds, PDF, CDF, CCDF and box plot[START_REF] Mcgill | Variations of box plots[END_REF][START_REF] Hosking | Parameter and quantile estimation for the generalized pareto distribution[END_REF][START_REF] Williamson | The box plot: a simple visual method to interpret data[END_REF][START_REF] Osterwald-Lenum | A note with quantiles of the asymptotic distribution of the maximum likelihood cointegration rank test statistics[END_REF][START_REF] Helton | Probability, conditional probability and complementary cumulative distribution functions in performance assessment for radioactive waste disposal[END_REF]. In general, the last four presentation patterns are usually preferable to the first several indices which will make large amount of uncertainty information neglected in implementation. Furthermore, box plot is definitely beneficial for displaying the uncertainty in model output with normative input factors and comparing the uncertainties in a number of related variables. The box plot is a standardised way of displaying the distribution of data. It is often generated by a box and whisker plots. The bottom and top of the box are, in general, the first and third quartiles of all of the data. The band inside the box is always the second quartile (the median).

Figure 1 .= x 2 1 + x 2 2 . ( 1 . 3 )Figure 1 . 3 :

 1221313 Figure 1.3: Representation of uncertainty in model output. (a): uncertainty analysis for the model presented in Eq. (1.3) whose input uncertainty is defined by Fig. 1.2(a). (b): box plot for the normalized susceptible and recovered agents at equilibrium state of SIR model with input factors assumed to be uniformly distributed between 0 and 1. Bars show the full range of the ensemble distribution of values; boxes show the range encompassed by the 25th and 75th percentiles; the horizontal line and square within each box show the median and mean, respectively.

Figure 1 . 4 :

 14 Figure 1.4: Progression of population for SIR model.

Figure 1 .Figure 1 . 5 :

 115 Figure 1.5: Representation of uncertainty in model output which acts a function of time parameter. Analysis of the time-dependent normalized number of infected agents generated in SIR model is presented as a particular example. Three involved input factors are assumed to be uniformly distributed between 0 and 1.

  ) To determine the optimum settings of input factors that result in the maximum or minimum output value over a certain range of interest. In general, the first-degree model is sufficient to determine which of input factors affect the model output of interest most. For a deep understanding of the effects produced by input factors on the model output, however, a more complicated design should be implemented to estimate a second-degree polynomial model. A series of experiments should be first designed to perform response surface analysis, helping generate the mapping from input factors to the output variable. The design, denoted by D, can be represented by an M × n matrix, as displayed in Eq. (1.1), where M is the number of experiments (the size of a design) and n the number of input variables. Each row of D represents a point in the n-dimensional input space. Designs used for estimating the first-degree model are usually referred to as first-order designs and those used for estimating the second-degree model, second-order designs.

1 Figure 1 . 6 :

 116 Figure 1.6: Example of random sampling to generate a sample of size M = 20 from x = (x1, x2) with x1 uniformly distributed and x2 normally distributed, defined by Fig.1.2(a). First two panels illustrate sampling of values for x1 and x2, respectively, and the last panel presents a random pairing of the sampled values os x1 and x2 in the construction of a random sampling.

Figure 1 .

 1 6 presents an example of random samples of size M = 20 generated from x = (x 1 , x 2 ). x 1 is uniformly distributed in the real range [0, 1] and x 2 follows the standard normal distribution (Fig.1.2(a)). Clustering effect and empty areas are evident in the distribution of sample points. The uniformity of sample points in input space directly determines the accuracy of analysis results when the number of sample points is finite.

1 Figure 1 . 7 :

 117 Figure 1.7: Example of LHS to generate a sample of size M = 20 from x = (x1, x2) with x1 uniformly distributed and x2 normally distributed, defined by Fig.1.2(a). First two panels illustrate sampling of values for x1 and x2, respectively, and the last one a random pairing of sampled values of x1 and x2 in the construction of a LHS.

Figure 1 .

 1 Figure 1.7 presents an example of LHS of size 20 for a special two-variate situation x = (x 1 , x 2 ). x 1 is uniformly distributed between 0 and 1 and x 2 follows the standard normal distribution, as defined in Fig.1.2(a).

1 Figure 1 . 8 :

 118 Figure 1.8: An example of QMC sampling based on Sobol sequences. The sample size: M = 20. x1 is uniformly distributed and x2 normally distributed, defined by Fig.1.2(a). The first two panels illustrate sampled values of x1 and x2, respectively, and last one the distribution of sampled values in the construction of QMC sampling.

  example of Sobol sequences-based QMC sampling is presented in Fig.1.8. A sample of size M = 20 is generated from x = (x 1 , x 2 ). x 1 is assumed to be uniformly distributed between 0 and 1, and x 2 normally distributed (Fig.1.2(a)).A qualitative comparison is also presented in Fig.1.9 for the uniformity of sampled values provided by different sampling strategies, including random sampling, standard LHS, LHS with the use of maximin distance criterion, LHS with the use of CL 2 criterion, and QMC sampling based on Sobol sequences. A situation of 2-dimensional input vector is exemplified. Two kinds of sample size: M = 16 and M = 100 are considered, respectively. Results suggest a higher uniformity of sampled values provided by Sobol sequences-based QMC sampling, compared to those given by random sampling and LHSs. Furthermore, the optimised LHS is superior to the standard LHS which precedes random sampling in constructing uniformly distributed sample values, especially when a large number of input factors are included. Some other newly introduced sampling

  labels the output vector. It may contain several independent or correlated observations. For simplicity in statement, we just consider an arbitrary element in y, denoted by y. The mapping from input vector to the model output can be analytically obtained if the functional form of the model under study is previously provided. With this kind of situation, optionally, the mapping of interest can also be constructed by computational strategies (sampling-based methods) which are most frequently used since their implementations do not require the form of model function to be given at first. In the following, several representations of analysis results are introduced by the use of sampling-based methods.1.4.1 Scatter plotsScatter plots-based representation is a natural starting point in the uncertainty and sensitivity analysis of complex systems. It provides an intuitive visual indication of the dependence of output variable upon individual input factors. Furthermore, it is definitely advantageous for directly understanding the impact of uncertainties in input factors on the uncertainty in model output. Scatter plots of two two-variate simple models are presented as examples: one is constructed in the absence of coupling items between input factors (Fig. 1.10) and the other is designed in the presence of coupling items. Input factors are assumed to be independent of each other and to be uniformly distributed in the real range [0, 1]. The scatter plots generated by the standard LHS and by the one-at-a-time method are considered, respectively. For sampling-based methods, the dependence of output y versus a single input factor is plotted in keeping values of the rest factors arbitrary. Regarding one-at-a-time method, however, the relationship of y (a) Sampled values based on random sampling. (b) Sampled values based on standard LHS. (c) Sampled values based on LHS with the use of maximin distance criterion. (d) Sampled values based on LHS with the use of CL 2 criterion. (e) Sampled values of QMC sampling based on Sobol sequences.

Figure 1 . 9 :

 19 Figure 1.9: Distribution of sampled values in a 2-dimensional input space. Input factors are assumed to be uniformly distributed in [0, 1]. (a): sampled values are generated by employing random sampling based on pseudo-random numbers. (b): sampled values are generated by applying the standard LHS. (c): sampled values are generated based on LHS which is optimised by the maximin distance criterion. (d): sample values are generated based on LHS which is optimised by CL2 criterion (e): sampled values are generated by Sobol sequences-based QMC sampling. From left to right, the size of sample M = 16 and M = 100, respectively. The unit square is divided into 64 (on the left) and 256 (on the right) squares.

Figure 1 . 10 :

 110 Figure 1.10: (Model 1) Examples of scatter plots-based representation for the uncertainty and sensitivity analysis of model: y = 4x1 + 2x 2 2 . The input factors are assumed to be independent of each other and to be uniformly distributed in the real range [0, 1]. Panel (a) presents the plots generated by standard LHS. Panel (b) illustrates the plots obtained in terms of one-at-a-time method.

Figure 1 . 11 :

 111 Figure 1.11: (Model 2) Examples of scatter plots-based representation for the uncertainty and sensitivity analysis of model: y = 4x1 + 2x 2 2 -4x1x2. The input factors are assumed to be independent of each other and to be uniformly distributed in the real range [0, 1]. Panel (a) presents the plots generated by standard LHS. Panel (b) illustrates the plots obtained in terms of one-at-a-time method.

Figure 1 .

 1 Figure 1.12: (SIR model) Examples of scatter plots-based representation for the uncertainty and sensitivity analysis of SIR model. The input factors β (infectious probability), γ (recovered rate), and s(0) (initial proportion of susceptible agents) are assumed to be independent of each other and uniformly distributed in the real range [0, 1]. Panel (a) presents the plots generated by standard LHS. Panel (b) illustrates the plots obtained in terms of one-at-a-time method.

(Fig. 1 .

 1 10(a)) (Fig. 1.11(a)) (Fig. 1.12(a)) (Fig. 1.10(b)) (Fig. 1.11(b)) (Fig. 1.12(b))

(Fig. 1 .

 1 10(a)) (Fig. 1.11(a)) (Fig. 1.12(a)) (Fig. 1.10(b)) (Fig. 1.11(b)) (Fig. 1.12(b))

(Fig. 1 .

 1 10(a)) (Fig. 1.11(a)) (Fig. 1.12(a)) (Fig. 1.10(b)) (Fig. 1.11(b)) (Fig. 1.12(b))

  It provides a global measure of the strength of the association between input variables and the output variable. The coefficient of multiple correlation takes values between 0 and 1. A higher value indicates a stronger linear association connecting input factors and output variable. A value of zero suggests the absence of linear relationship between input factors and output variable. Unlike the previously introduced correlation coefficients which offer both the strength and direction of the underlying association connecting input factors and the output varaible, the coefficient of multiple correlation tells only the strength of the association. Regarding a general model of two input factors, the coefficient of multiple correlation is defined by

  Amatrices is presented in table 1.6 with sample size M = 5 and three input variables considered. They are assumed to be independently and uniformly distributed in the real range [0, 1].

  x 0 and x m denote the lower and upper bounds of x, respectively. Detailed derivation of M k (x) is presented in Appendix A. Inserting the expression of central moments (Eq. (2.13)) into Eq.(2.8)

Figure 2 . 1 :

 21 Figure 2.1: Distributions of the observation R γ v with different parameters for power-law function: y = x α (panel (a)) and exponential function: y = e αx (panel (b)) by assuming input variable is uniformly distributed; To the left in both panels is the distribution of R γv versus the order index γ for two kinds of nonlinear relationship connecting y and x with constant values of parameter α; to the middle is the dependence of R γ v on the parameter α in considering different order contributions of uncertainty in x; to the right is the association of R γ v with the distribution parameter of variable x by setting α = 4 for both panels. For a function with specified parameter α, the observation R γ v only depends upon the ratio of µ to σ (power-law function) or σ (exponential function) for each specified γ. In the left and middle plots, µ/σ = 2 in panel (a) and σ 2 = 0.5 in panel (b).

Figure 2 . 2 :

 22 Figure 2.2: Distributions of the observation R γ v versus different parameters for power-law function: y = x α (panel (a)) and exponential function: y = e αx (panel (b)) by assuming input variable is normally distributed; To the left in both panels is the distribution of R γ v with order index γ for two kinds of functional form connecting y and x, by considering different values of parameter α; to the middle is the association of R γ v with parameter α under considering different order contributions of uncertainty in x; to the right is the dependence of R γ v upon the distribution parameter of variable x by setting α = 5 in panel (a) and α = 1 in panel (b), respectively. For a function with specified α, the observation R γ v correlates with the ratio of µ to σ (power-law function) or with σ (exponential function) for each fixed γ. In the left and middle plots, µ/σ = 2 in panel (a) and σ = 1 in panel (b).

  .2. Panel (a) illustrates the data for the power-law relationship connecting y and x. Panel (b)

Figure 2 . 3 :

 23 Figure 2.3: Results associated with the observation R γ v for exponential relationship connecting y and x: y = e αx by assuming input x is uniformly and normally distributed, respectively. To the left is the distribution of R γ v versus order index γ by keeping α = σ = 1; To the middle denotes the dependence of R γ v upon parameter α in fixing γ = σ = 1; To the right is the correlation between R γ v and σ by setting γ = α = 1.

. 35 )Figure 2 . 4 :Figure 2 . 5 :

 352425 Figure 2.4: Typical wind turbine power output with steady wind speed.

Figure 2 . 6 :

 26 Figure 2.6: A comparison of our method with Sobol's one. Sobol's values are obtained by random sampling-based strategy (sample size: M = 10000), with error bars indicating the standard deviation of 100 runs.

Figure 2 . 7 :

 27 Figure 2.7: A comparison of our method with Sobol's one. Sobol's values are obtained by random sampling-based strategy (sample size: M = 10000), with error bars indicating the standard deviation of 100 runs.

  C -1 is the inverse of matrix C, µ(r) = (µ(r 1 ), µ(r 2 ), • • • , µ(r n )) T the mean vector of newly introduced random variables, µ = (µ 1 , µ 2 , • • • , µ n ) T the mean vector of input variables, and similarly for σ(r) and σ the vectors of standard deviations. Entries in matrices A and C: {a ij , c ii ; i, j = 1, 2, • • • , n} are called coefficients that specify the correlated and independent sections of input variables. They are determined by given pairwise correlations through equations

40 )

 40 By the use of Eq. (3.4), we get the exact expression of the variance of model response:V (y) = 9 + 4ρ 12 + 8ρ 13 + 4ρ 23 ,(3.41)

Figure 3 . 1 :

 31 Figure 3.1: The convergence of analysis results with the order of Taylor expansion.

  3 lists the exact analytic values of the underlying sensitivity indices. Data show a vanishing coupling effect between input independence and correlations. This because the nonlinear problems of single variables are absent in the form of model function. The convergence of our analysis results is displayed in Fig. 3.2, along the direction of the highest-order of Taylor expansion that is considered in our method. Values are convergent at γ = 2 (the highest-order of Taylor expansion of the model output).

Figure 3 . 2 :

 32 Figure 3.2: The convergence of analysis results with the order of Taylor expansion.

Figure 3 . 3 :

 33 Figure 3.3: The convergence of analysis results with the order of Taylor expansion.

(

  indicated by the standard deviation) of R 0 for independent and correlated situations, respectively. A ranking of input parameters is displayed in Fig.3.4, according to the total sensitivity indices. Two values of γ (the highest-order of Taylor expansion that is considered in our method) are considered. Results suggest that the analysis results of γ = 1 are almost the same as those of γ = 2 for both the absence and presence of input correlations. This implies our analysis results are convergent at the first-order Taylor expansion for the model under discussion. In both the absence and presence of input correlations, parameters κ (rate of development to AIDS), β 2 (contact rate of susceptible with symptomatic infective), and n 2 (number of sexual partners of susceptible with symptomatic infective) are identified as the most influential, while parameters β 0 (the birth rate of infective) and γ (the fraction of susceptible newborn from infective class) are identified as negligible, in determining the basic reproduction number R 0 . This provides one an opportunity to effectively limit the spread of a disease by controlling the three most influential parameters, and to simplify an HIV model by neglecting the effects contributed by parameters β 0 and γ.

Figure 3 . 4 :

 34 Figure 3.4: (Color online) Parameters are ranked according to their global importance (the total sensitivity indices sT i) in establishing the uncertainty of the basic reproduction number R0 for HIV-1 model with vertical transmission. The left panel is for independent situation, and the right one for the correlated one with ρ β 1 n 1 = 0.3 and ρ β 2 n 2 = 0.5.

ρ β 1 n 1 =0Figure 3 . 5 :

 135 Figure 3.5:The distribution of our analysis results along the direction of the correlation between β1 and n1, and the direction of the correlation between β2 and n2. The sensitivity indices for n1 and n2 are exactly the same as those for β1 and β2, respectively.

  tributed in the closed interval [0, 1] ([0.1, 1.1] for the last case). The underlying results obtained by our analytic method and by the sampling-based Sobol's one are exhibited in tables 4.1 and 4.2, respectively. Two sampling strategies are employed here, including LHS and QMC sampling. Typically, QMC sampling is performed by using high-dimensional Sobol low-discrepancy sequences that hold advanced uniformity properties in the input space. In practice, LHS-based results are averaged over 1000 configurations used for eliminating the randomness in sample generation. And QMC sampling-based values are obtained by unique configuration as low-discrepancy sequences are deterministic for given sample size and input variables. Results indicate a good agreement between our values and Sobol's ones, except the last two examples involving interaction effect of nonlinear problems of x 3 . For these two cases, Sobol's sensitivity indices overestimate the impacts contributed by x 1 , x 2 , as well as by their interaction, but underestimate the roles played by the association of each of them with x 3 , compared to the values provided by our method. This difference is introduced by the following components separately appearing in the analytic expansions of s 13 , s 23 , and s 123 (see the definition of sensitivity indices in the second chapter):

  This was soon followed by a deterministic model, i.e., SIR model, proposed by W.O. Kermack and A.G. McKendrick in 1927 with the consideration of a fixed population divided into three compartments: susceptible group S who are not infected with disease but can get infected, infectious category I who are infected with disease and can spread the disease to group S, and recovered community R who have been recovered from infection and are immune to reinfection [37]. Another widely considered model for epidemic spreading is SIS model relying on a coarse grained description of individuals in the population within which individuals are immediately susceptible once they have recovered.

(4. 15 )

 15 Numerical simulations of Eqs (4.9) and (4.14) are presented in Fig. 4.1 by considering different values of input variables. Results demonstrate that the underlying system will fast reach the equilibrium state for both models. Uncertainty and sensitivity analysis of both models is performed by sampling-based strategies including LHS and QMC sampling. QMC sampling is implemented by the use of high-dimensional Sobol lowdiscrepancy sequence generator. The underlying analysis results are exhibited in table 4.3 where rows of SIR r and SIS r correspond to the values after introducing our proposed modification. LHS-based values are averaged over 1000 configurations for eliminating the randomness in sample generation, and QMC sampling-based ones are provided by unique configuration.

(4. 20 )

 20 By introducing modifications presented above, analysis results of SIR model suggest that the initial proportion s(0) of susceptible individuals plays the most important role in the estimation of equilibrium state of epidemic spreading; infectious probability β and recovered rate γ almost explain the same proportion of variance in response s; the interaction effects of input factors play fragile roles in the decision of model response.Regarding SIS model, parameter γ plays more robust role than β, and their interaction effect also provides a crucial contribution to the uncertainty of equilibrium state, by holding 31% of the decision of s. Sampling-based plots for the association of s with each single factor are exhibited in Fig.4.3, together with phase diagrams of model response s in any possible two-dimensional input spaces.

  Numerical simulation for SIS model

Figure 4 . 1 :

 41 Figure 4.1: Numerical simulation results for both SIR and SIS models by considering different values of input parameters.

Figure 4 . 2 :

 42 Figure 4.2: The normalized number s of susceptible individuals acts as functions of individual input factors at equilibrium state of SIR model.

Figure 4 . 3 :

 43 Figure 4.3: Scatter plots based on LHS, as an example, with M = 1000 for the normalized number of susceptible individuals at equilibrium state, together with its phase diagrams in any possible two-dimensional input space.

Figure 5 . 1 :

 51 Figure 5.1: Time evolution of opinion dynamics for different specified values of parameters with N = 1000 and ka = 4. One unit time corresponds to sampling one round of gambling game.

Figure 5 . 2 :

 52 Figure 5.2: The generation of three adjacent regions in the opinion space.

5. 4 . 1

 41 Theoretical analysis for special cases p = 0, ε = 2

Fig. 5 Figure 5 . 3 :

 553 Fig. 5.3(c).

Figure 5 . 4 :

 54 Figure 5.4: Show only single run with N = 1000, ka = 4. Model parameters: p = 0, ε = 2, γ = 0.003. One unit time corresponds to sampling one round of gambling game.

Figure 5 . 5 :

 55 Figure 5.5: The generation of four adjacent regions in the opinion space.

Figure 5 . 6 :

 56 Figure 5.6: A comparison between numerical results (star symbols) and analytic values (solid line) for the association of Ωs with parameter ε in keeping p = γ = 0. Numerical results are averaged over 1000 simulations with N = 1000 and ka indicated on the figure. Analytic results are represented, from left to right, by Eqs. (5.39), (5.40) and (5.41).

Figure 5 . 7 :

 57 Figure 5.7: Steady values of order parameter Ωs versus confidence threshold ε with the impacts of mindless agents in keeping resource undivided (γ = 0) (a) and the roles of resource allocation parameter in keeping all agents economic (p = 0) (b) considered, respectively. In panel (a), the fraction p of mindless agents is 0, 0.4, 0.6 and 0.8 in sequence from up to down. In panel (b), the resource allocation parameter γ = 0, 0.1, 0.3, 0.5 and 0.7 from up to down.Results are averaged over 1000 simulations with N = 1000 and ka = 4.

  Fig. 5.7.

Figure 5 . 8 :

 58 Figure 5.8: Show only single runs. Opinion evolution with time step t for N = 100, ka = 4. The initial state of opinion distribution keeps the same for four considered cases. One unit time corresponds to one sampling round of gambling game.

. 45 )

 45 This critical value of confidence threshold indicates the emergence of consensus state in keeping p = γ = 0.When parameter γ is limited to the range (0, 1/3], initial winners must hold opinions no larger than -(1 + γ)/4. By changing 1/4 to (1 + γ)/4 in Eq. (5.28), Eq. (5.29) is then updated by

  ) where M labels the number of clusters, m j the normalized size of jth cluster. That Y → 0 as N → ∞ indicates random opinion distribution in the range [-1, 1] since Y = 1/M labels M clusters of equal size. Finite Y corresponds to the emergence of large clusters. The generation of large cluster can be identified by the relative difference between quantity Y and the largest cluster, given by R(Y, m L ) = Y -(m L ) 2 Y . (5.53) R(Y, m L ) ≈ 0 indicates Y just considers the contribution of m L compared to which other clusters, if exist, can be neglected, and R(Y, m L ) ≈ 1, which is equivalent to Y ≈ 0, indicates a vanishing large cluster but existent small ones with comparable size. The underlying results regarding Y and R(Y, m L ) are exhibited along ε-axis (Fig. 5.9).

Figure 5 . 9 :Figure 5 . 10 :

 59510 Figure 5.9: Steady values of dispersion index Y and relative difference R(Y, m L ) versus parameter ε, for different values of p and γ indicated on the figure. Network parameters used for the results generation averaged over 1000 simulations are N = 1000 and ka = 4.

  Phase diagrams of n w generated in γ-ε space for p=0.2, 0.5, and 0.8 from left to right. n e w as a function of (1-p) with γ = 0 to the left and ε = 2 to the right.

Figure 5 . 11 :

 511 Figure 5.11: (Color online) Distribution of the fraction nw of winners (a), and the dependence of the percentage of economic agents in winners upon the percentage of economic agents included in the system under consideration (b). Data are averaged over 100 simulations for panel (a) and 1000 simulations for panel (b) with N = 1000 and ka = 4.

23 r 1 .

 1

V

  [START_REF] Lee | Response surface augmented moment method for efficient reliability analysis[END_REF] can be equivalently obtained by formulating x 3 on the bases of x 1 and x 2 , constituted of

  In the second nonlinear test case presented in Chapter 3, the fourth-order central moment with four different variables is involved:M 1,1,1,1 (x 1 , x 2 , x 3 , x 4 ) = E [(x 1 -µ 1 )(x 2 -µ 2 )(x 3 -µ 3 )(x 4 -µ 4 )] ,(D.1)

(E. 6 )V 2 7 )

 627 By inserting Eqs. (A.5) and (E.6) into Eq. (E.5), we get Zero mean value of x 1 also suggests a vanishing main partial variance contribution produced by x 3 alone since x 3 just appears in the combination with sin function of x 1

r 1 = ρ 13 x 3 + √ 1 -ρ 2 13 r 1 ,(E. 8 ) with r 1

 11181 holding the same distribution as x 1 (since the mean values of x 1 and x 3 are zero). Equation (E.4) is then rewritten as

-ρ 2 13 )ρ 2 13 )ρ 2 13 )ρ 2 13 )ρ 2k 13 ( 1 -ρ 2 13 )

 1313131313113 i+j+1-k (2i + 2j + 3 -2k)(2k + 1). (E.10)Firstly, consider the generation of x 1 on the basis of x 3 . Inserting Eq. (E.8) into covariance items in Eq. (E.14) yields i+j+1-k M 2i+2j+2-2k (r 1 )M 8+2k (x 3 ) i+j+1-k (2(i + j + 1) -2k + 1)(8 + 2k + 1) , i+j+1-k M 2(i+j+1)-2k (x 1 )M 4+2k (x 3 ) i+j+1-k(2(i + j + 1) -2k + 1)(4 + 2k + 1)

V 13 =(- 1 )

 131 i+j π 2(i+j+1)+4 (2i + 1)!(2j + 1)! • ( 2(i + j + 1)2k)

r 3 = ρ 13 x 1 + √ 1 -ρ 2 13 r 3 ,r 3

 3133 holding the same distribution as x 3 . Then we get -ρ 2 13 ) 4-k M 2(i+j+1)+2k (x 1 )M 8-2k (x 3 ) ρ 2 13 ) 2-l M 2(i+j+1)+2l (x 1 )M 4-2l (x 3 ) i+j+1)+4 ρ 2l 13 (1 -ρ 2 13 ) 2-l (2(i + j + 1) + 2l + 1)(4 -2l + 1), 2j + 2l + 3)(5 -2l) ]

Figure F. 1 :

 1 Figure F.1: The generation of four adjacent regions in the opinion space when ε ∈ [0, 1/4].

Figure F. 2 :

 2 Figure F.2: The generation of four adjacent regions in the opinion space when ε ∈ (1/4, 1/2].

Figure F. 3 :

 3 Figure F.3: The generation of five adjacent regions in the opinion space when ε ∈ (1/2, 3/4].

Figure F. 4 :

 4 Figure F.4: The generation of three adjacent regions in the opinion space when ε ∈ (5/4, 7/4].

29 )

 29 Consider sufficiently high confidence threshold: ε > 7/4. The partition of opinion space can be simplified as two regions (Fig. F.5). With local communication, almost all agents

Figure F. 5 :

 5 Figure F.5: The generation of two adjacent regions in the opinion space when ε ∈ (7/4, 2].

Table 1 . 1 :

 11 Estimated mean value and its standard deviation based on MC and importance sampling (denoted by IS). Results are obtained by N independent simulations of sample size 100.

	N	Simulation	Ê(f ; P ) σ( Ê(f ); P )
	100	MC IS	0.459 0.460	0.051 0.020
	1000	MC IS	0.459 0.461	0.048 0.021
	100000	MC IS	0.460 0.460	0.049 0.022

Table 1 . 2 :

 12 The CCs for two-variate nonlinear models and SIR model with scatter plots presented in Figs. 1.10, 1.11 and 1.12. For SIR model, x1 denotes the parameter β, x2 the parameter γ, x3 the parameter s(0), and y the output variable s.

		sampling-based method	one-at-a-time method
	CC	model 1	model 2	SIR model	model 1	model 2	SIR model

.68) As an example, the CCs are shown in table 1.2 for two two-variate nonlinear models and SIR model. Their scatter plots are presented in Figs. 1.10, 1.11 and 1.12, by the use of LHS-based method and one-at-a-time method.

The partial correlation coefficient (PCC) provides a measure to characterise the degree of association between a single factor x i and an output variable y with the effects of the remaining elements of x removed. By the use of linear regression model, a new variable

Table 1 . 3 :

 13 RCCs for two-variate nonlinear models and SIR model with scatter plots presented in Figs. 1.10, 1.11 and 1.12. In SIR model, x1 denotes the parameter β, x2 the parameter γ, x3 the parameter s(0), and y the output variable s.

	RCC	sampling-based method model 1 model 2 SIR model	one-at-a-time method model 1 model 2 SIR model

Table 1 . 4 :

 14 KRCCs for two-variate nonlinear models and SIR model with scatter plots presented in Figs. 1.10, 1.11 and 1.12. For SIR model, x1 denotes the parameter β, x2 the parameter γ, x3 the parameter s(0), and y the output variable s.

	KRCC	sampling-based method model 1 model 2 SIR model	one-at-a-time method model 1 model 2 SIR model

Table 1 . 5 :

 15 The PCCs and PRCCs for two-variate nonlinear models and SIR model with scatter plots presented in Figs. 1.10(a) (model 1), 1.11(a) (model 2) and 1.12(a) (SIR model). For SIR model, x1 denotes the parameter β, x2 the parameter γ, x3 the parameter s(0), and y the output variable s.

	τ (x1, y)	0.71	0.73	-0.27	1	1	-1
	τ (x2, y)	0.29	0.05	0.46	1	0.00	1
	τ (x3, y)	-		-	0.45	-	-	1
	measures	PCC model 1 model 2 SIR model model 1 model 2 SIR model PRCC
	ρ(x1, y)		0.99	0.83	-0.57	0.98	0.90	-0.68
	ρ(x2, y)		0.97	-0.02	0.64	0.92	0.17	0.78
	ρ(x3, y)		-	-	0.73	-	-	0.76

As an example, the KRCCs are exhibited in table 1.4 for the models with scatter plots presented in Figs. 1.10, 1.11 and 1.12, by the use of LHS-based method and one-at-atime method. Similar to PCC, the partial rank correlation coefficient (PRCC) performs a partial correlation on rank-transformed data: rx i and ry by the use of linear regression models described by Eqs. (1.69) and (1.70). The PRCC provides a robust sensitivity measure for nonlinear but monotonic relationships between x i and y in the absence of input correlations. The PCCs and PRCCs are presented in table 1.5 for nonlinear models with scatter plots illustrated in Figs. 1.10(a) (model 1), 1.11(a) (model 2) and 1.12 (SIR model), by the use of LHS-based method.

Table 1 . 6 :

 16 An example of construction of both AA matrices with sample size M = 5 and number of input variables n = 3. Matrices A and B are generated by the use of standard LHS. Three input variables are independently and uniformly distributed in the range [0, 1].

	(i) B and B	(i)

Table 1 . 7 :

 17 Variance-based sensitivity analysis results for the model of form y = 4x1 + 2x 2 2 . The standard LHS and QMC sampling are used, with the effect of sample size M considered. Both Sobol's (Eqs. (1.80)-(1.81)) and Jansen's (Eqs. (1.85)-(1.86)) methods are employed to the evaluation of sensitivity indices. For LHS, analysis results are averaged over 100 independent configurations, with standard deviations displayed in parentheses. QMC sampling, generated based on Sobol sequences, is deterministic and numbers in parentheses indicate the thousandth and ten-thousandth digits.

	sampling	M	method	E(y)	V (y)	s 1	s 2	s T 1	s T 2
		10 2	Sobol' Jansen	2.67(.00) 2.67(.00)	1.69(.15) 1.70(.15)	0.79(.05) 0.79(.03)	0.21(.08) 0.22(.11)	0.79(.08) 0.78(.11)	0.21(.05) 0.21(.03)
	LHS
		10 4	Sobol' Jansen	2.67(.00) 2.67(.00)	1.69(.01) 1.69(.01)	0.79(.00) 0.79(.00)	0.21(.01) 0.21(.01)	0.79(.01) 0.79(.01)	0.21(.00) 0.21(.00)
		10 2	Sobol' Jansen	2.66(72) 2.66(72)	1.67(50) 1.67(50)	0.79(63) 0.77(61)	0.21(41) 0.22(85)	0.78(60) 0.77(15)	0.20(37) 0.22(40)
	QMC	10 4	Sobol' Jansen	2.66(67) 2.66(67)	1.68(89) 1.68(89)	0.78(90) 0.78(95)	0.21(05) 0.21(06)	0.78(95) 0.78(95)	0.21(10) 0.21(05)
		3 • 10 4	Sobol' Jansen	2.66(66) 2.66(66)	1.68(88) 1.68(88)	0.78(97) 0.78(95)	0.21(03) 0.21(05)	0.78(97) 0.78(95)	0.21(04) 0.21(05)

Table 1 . 8 :

 18 Variance-based sensitivity analysis results for the model of form y = 4x1 + 2x 2 2 -4x1x2. The standard LHS and QMC sampling are used, with the effect of sample size M considered. Both Sobol's (Eqs. (1.80)-(1.81)) and Jansen's (Eqs. (1.85)-(1.86)) methods are employed to the evaluation of sensitivity indices. For LHS, analysis results are averaged over 100 independent configurations, with standard deviations displayed in parentheses. QMC sampling, generated based on Sobol sequences, is deterministic and numbers in parentheses indicate the thousandth and ten-thousandth digits.

	sampling	M	method	E(y)	V (y)	s 1	s 2	s T 1	s T 2
		10 2	Sobol' Jansen	1.67(.04) 1.67(.03)	0.46(.04) 0.47(.04)	0.71(.19) 0.72(.05)	0.03(.21) 0.04(.14)	0.97(.21) 0.96(.12)	0.29(.19) 0.29(.04)
	LHS
		10 4	Sobol' Jansen	1.67(.00) 1.67(.00)	0.47(.00) 0.47(.00)	0.72(.02) 0.72(.00)	0.05(.02) 0.05(.02)	0.95(.02) 0.95(.01)	0.28(.02) 0.28(.00)
		10 2	Sobol' Jansen	1.67(46) 1.67(46)	0.48(74) 0.48(74)	0.67(71) 0.72(51)	0.08(13) 0.09(33)	0.91(87) 0.85(59)	0.32(29) 0.26(22)
	QMC	10 4	Sobol' Jansen	1.66(67) 1.66(67)	0.46(65) 0.46(65)	0.71(46) 0.71(41)	0.04(78) 0.04(73)	0.95(22) 0.95(24)	0.28(54) 0.28(60)
		3 • 10 4	Sobol' Jansen	1.66(67) 1.66(67)	0.46(68) 0.46(68)	0.71(41) 0.71(43)	0.04(72) 0.04(77)	0.95(28) 0.95(25)	0.28(59) 0.28(57)

Table 1 .

 1

	sampling	M	method	E(y)	V (y)	s β	sγ	s s(0)	s T β	s T γ	s T s(0)
		10 2	Sobol' Jansen	0.24(.01) 0.25(.01)	0.07(.01) 0.07(.01)	0.19(.14) 0.19(.12)	0.20(.12) 0.19(.13)	0.43(.15) 0.42(.10)	0.30(.09) 0.31(.07)	0.33(.12) 0.33(.08)	0.59(.13) 0.60(.09)
	LHS
		10 4	Sobol' Jansen	0.24(.00) 0.24(.00)	0.07(.00) 0.07(.00)	0.19(.01) 0.19(.01)	0.18(.01) 0.19(.01)	0.42(.01) 0.42(.01)	0.31(.01) 0.31(.01)	0.34(.01) 0.34(.01)	0.59(.01) 0.59(.01)
		10 2	Sobol' Jansen	0.23(83) 0.23(83)	0.07(30) 0.07(30)	0.23(21) 0.31(77)	0.23(95) 0.36(61)	0.31(50) 0.41(09)	0.34(66) 0.32(83)	0.50(42) 0.41(59)	0.54(42) 0.56(16)
	QMC	10 4	Sobol' Jansen	0.24(19) 0.24(19)	0.07(14) 0.07(14)	0.18(56) 0.18(30)	0.18(11) 0.18(08)	0.41(54) 0.41(36)	0.30(97) 0.30(99)	0.33(99) 0.33(97)	0.59(45) 0.59(50)
		3 • 10 4	Sobol' Jansen	0.24(19) 0.24(19)	0.07(14) 0.07(14)	0.18(61) 0.18(75)	0.18(15) 0.18(30)	0.41(51) 0.41(59)	0.30(82) 0.30(81)	0.34(13) 0.34(08)	0.59(53) 0.59(55)

9: Variance-based sensitivity analysis results for SIR model. The standard LHS and QMC sampling are used, with the effect of sample size M considered. Both Sobol's (Eqs. (1.80)-(1.81)) and Jansen's (Eqs. (1.85)-(1.86)) methods are employed to the evaluation of sensitivity indices. For LHS, analysis results are averaged over 100 independent configurations, with standard deviations displayed in parentheses. QMC sampling, generated based on Sobol sequences, is deterministic and numbers in parentheses indicate the thousandth and ten-thousandth digits.

Table 1 .

 1 10: An example of the construction of two independent QMC samplings with sample size M = 10 and number of input variables n = 3. Three input variables are independently and uniformly distributed in the real range [0, 1].

	              	0.09375 0.59375 0.84375 0.34375 0.46875 0.96875 0.71875 0.21875 0.15625 0.65625	0.46875 0.96875 0.21875 0.71875 0.09375 0.59375 0.34375 0.84375 0.15625 0.65625	0.84375 0.34375 0.09375 0.59375 0.46875 0.96875 0.71875 0.21875 0.53125 0.03125	0.40625 0.90625 0.15625 0.65625 0.28125 0.78125 0.03125 0.53125 0.84375 0.34375	0.28125 0.78125 0.53125 0.03125 0.65625 0.15625 0.40625 0.90625 0.84375 0.34375	0.34375 0.84375 0.09375 0.59375 0.71875 0.21875 0.96875 0.46875 0.40625 0.90625

Table 2 . 1 :

 21 Parameters associated with the uncertainty in wind speed x.

		µx	V (x) M3(x) M4(x)
	c = 4, k = 2	4.79 1.97	3.03	16.09
	c = 4, k = 4	4.1	0.5	0.21	0.73
	c = 10, k = 2 7.47 6.03	0.37	69.63

Table 2 . 2 :

 22 Analysis results of the power grid system in considering different values of γ (the highest order of Taylor expansion considered in our method).

				γ V (P d )/P 2 r	sx	sε	sxε
	c = 4, k = 2		1 0.01(35) 2 0.01(87)	0.26(07) 0.73(93) 0 0.46(51) 0.53(49) 0
	c = 4, k = 4		1 0.01(03) 2 0.01(05)	0.03(01) 0.96(99) 0 0.04(93) 0.95(07) 0
	c = 10, k = 2	1 0.08(28) 2 0.08(91)	0.87(92) 0.12(08) 0 0.88(77) 0.11(23) 0
	Sensitivity Index	0.5 1.0	c=4, k=2	c=10, k=2	c=4, k=4
		0.0
			x		x	x
						Sobol

Table 2 . 3 :

 23 Analysis results for the EOQ model while considering different values of γ.

	γ = V (Q * )	sD	sK	s h	sDK	s Dh	s Kh	s DKh
	1	2119	0.385 0.307 0.307	0	0	0	0
	2	2192	0.377 0.300 0.314 -0.002 0.006 0.005	0
	Input variables are assumed to be independently and uniformly distributed in real ranges
	as follows [225]:
		900 ≤ D ≤ 1600 unit per year,
		$75 ≤ K ≤ $125 per order,			(2.40)
		$1.5 ≤ h ≤ $2.5 per order and per year,
	which yields
	µ D = 1250,	µ K = 100,		µ h = 2,
	V (D) = 40833.333,	V (K) = 208.333,		V (h) = 0.083.	(2.41)
	tivity analysis are determined as
	s D = 0.377,	s K = 0.300,	s h = 0.314,
	s DK = -0.002,	s Dh = 0.006,	s Kh = 0.005,	s DKh = 0.000.	(2.42)

Substituting the functional form and distribution laws of model inputs into Eq.

(2.11)

gives the exact value of output uncertainty: V (Q * ) = 2195. The exact values of sensi-Sensitivity analysis results are listed in table 2.3 for different values of γ. While γ = 2,

  • • • ) is defined by Eq. (2.25), {µ} the central point in input space, and M i 1 ,••• ,in the multivariate central moment with n different variables. Mathematically, the multivariate central moment is defined as

Table 3 . 1 :

 31 Uncertainty and sensitivity analysis results for linear additive model by assuming different correlations between input variables.

	ρ	V (y)	x	si	s U i	s C i	s UC i	sT i	s TU i	s TC i	s TUC i
			x1 0.444 0.444	0.0	0.0	0.444	0.444	0.0	0.0
	ρ = 0	9	x2 0.111 0.111	0.0	0.0	0.111	0.111	0.0	0.0
			x3 0.444 0.444	0.0	0.0	0.444	0.444	0.0	0.0
			x1 0.328 0.118 0.210	0.0	0.590	0.118	0.472	0.0
	ρ12 = 0.8	12.2	x2 0.082 0.030 0.052	0.0	0.344	0.030	0.314	0.0
			x3 0.328 0.328	0.0	0.0	0.328	0.328	0.0	0.0
			x1 0.690 0.248 0.442	0.0	0.138	0.248 -0.110	0.0
	ρ12 = -0.8	5.8	x2 0.172 0.062 0.110	0.0	-0.379 0.062 -0.441	0.0
			x3 0.690 0.690	0.0	0.0	0.690	0.690	0.0	0.0
	ρ12 = 0.8		x1 0.225 0.072 0.152	0.0	0.629	0.072	0.557	0.0
	ρ13 = 0.5	17.8	x2 0.056 0.020 0.036	0.0	0.326	0.020	0.306	0.0
	ρ23 = 0.4		x3 0.225 0.169 0.056	0.0	0.539	0.169	0.371	0.0

where i ∈ {1, 2, 3}, j ∈ {1, 2}, k ∈ {1, 3}, l ∈ {2, 3}, and c ii , specifying the independence of variable x i , are determined via Eq. (3.17) as

Table 3 . 2 :

 32 Exact analytic results for uncertainty and sensitivity analysis of the first nonlinear model with different input correlations.

	ρ	V (y)	x	si	s U i	s C i	s UC i	sT i	s TU i	s TC

  ρ 12 ρ 13 ρ 14 ρ 12 1 ρ 23 ρ 24 ρ 13 ρ 23 1 ρ 34 ρ 14 ρ 24 ρ 34 1 = 12+4ρ 12 +4ρ 13 +8ρ 14 +2ρ 23 +4ρ 24 +4ρ 34 +ρ 2 13 +ρ 2 24 +2(ρ 12 ρ 34 +ρ 14 ρ 23 ), (3.55)

	
	
	
	 	.	(3.54)
	
	
	The total variance of model response is obtained by employing Eq. (3.4) as
	V (y)

Table 3 . 3 :

 33 Exact analytic values of uncertainty and sensitivity analysis for the second nonlinear model by assuming uncorrelated and correlated inputs.

	ρ	V (y)	x	si		s U i	s C i		s UC i	sT i		s TU i	s TC i	s TUC i
				x1 0.333 0.333	0.0		0.0	0.417 0.417	0.0	0.0
	ρ = 0	12		x2 0.083 0.083	0.0		0.0	0.167 0.167	0.0	0.0
				x3 0.083 0.083	0.0		0.0	0.167 0.167	0.0	0.0
				x4 0.333 0.333	0.0		0.0	0.417 0.417	0.0	0.0
	ρ13 = 0.5			x1 0.221 0.166 0.055		0.0	0.401 0.207	0.194	0.0
	ρ24 = 0.8	18.09	x2 0.055 0.020 0.035		0.0	0.323 0.040	0.283	0.0
				x3 0.055 0.041 0.014		0.0	0.235 0.083	0.152	0.0
				x4 0.221 0.080 0.141		0.0	0.489 0.100	0.389	0.0
	ρ12 = -0.5			x1 0.251 0.058 0.193		0.0	0.561 0.072	0.489	0.0
	ρ13 = 0.6	15.96	x2 0.063 0.030 0.033		0.0	0.0		0.060 -0.060	0.0
	ρ14 = 0.4			x3 0.063 0.024 0.039		0.0	0.298 0.049	0.249	0.0
				x4 0.251 0.148 0.103		0.0	0.514 0.185	0.329	0.0
	ρ12 = -0.5, ρ13 = -0.4		x1 0.289 0.110 0.179		0.0	0.208 0.137	0.071	0.0
	ρ14 = 0.2, ρ23 = 0.3	13.84	x2 0.072 0.032 0.040		0.0	0.150 0.065	0.085	0.0
	ρ24 = 0.4, ρ34 = 0.4			x3 0.072 0.040 0.032		0.0	0.179 0.080	0.099	0.0
				x4 0.289 0.121 0.168		0.0	0.699 0.152	0.547	0.0
		0.50					0.50
		0.25					0.25
	Sensitivity index	0.00 0.8					0.00 0.8
		0.4					0.4
		0.0					0.0
		0	2	4	6	8	10	0	2	4	6	8	10

Table 3 . 4 :

 34 Exact analytic values of uncertainty and sensitivity analysis for Ishigami function by assuming uncorrelated and correlated inputs. In case 1 x1 is generated based on x3, contrary to the second case where x3 is generated based on x1.

	ρ	V (y)	x	si	s U i	s C i	s UC i	sT i	s TU i	s TC i	s TUC i
			x1 0.036 0.036	0.0	0.0	0.557 0.557	0.0	0.0
	ρ13 = 0	13.845	x2 0.442 0.442	0.0	0.0	0.442 0.442	0.0	0.0
			x3	0.0	0.0	0.0	0.0	0.521 0.521	0.0	0.0
	ρ13 = 0.5		x1 0.039 0.044 0.039 -0.044 0.528 0.633 1.048 -1.153
	(case 1)	12.971	x2 0.472 0.472	0.0	0.0	0.472 0.472	0.0	0.0
			x3	0.0	-	-	-	0.489	-	-	-
	ρ13 = 0.5		x1 0.026	-	-	-	0.679	-	-	-
	(case 2)	19.110	x2 0.321 0.321	0.0	0.0	0.321 0.321	0.0	0.0
			x3	0.0	0.0	0.0	0.0	0.653 0.145 0.004	0.505
	model response owing to zero partial variance contributions associated with the inter-

action between x 2 and the rest. Consequently, we just consider here the correlation between x 1 and x 3 . The results of analytic analysis are listed in table 3.4 by assuming independent and correlated input variables. Two formulating strategies are considered in the presence of correlation: x 1 is formulated on the basis of x 3 and vice versa. They are non-equivalent for the uncertainty and sensitivity analysis of the model under discussion as x 1 and x 3 are uniformly distributed.

Table 3 . 5 :

 35 Description and baseline values of parameters for HIV/AIDS model, see Refs.[START_REF] Usman | Sensitivity analysis of the parameters of a model of hiv-1 transmission[END_REF][START_REF] Mukandavire | Hiv/aids model for assessing the effects of prophylactic sterilizing vaccines, condoms and treatment with amelioration[END_REF][START_REF] Safiel | Modelling the effect of screening and treatment on transmission of hiv/aids infection in a population[END_REF].

	Parameter	Symbol Baseline value
	Recruitment rate	Q0	0.029
	Birth rate of infective	β0	0.03
	Fraction of susceptible newborn from infective class	γ	0.4
	Contact rate of susceptible with asymptomatic infective	β1	0.2
	Contact rate of susceptible with symptomatic infective	β2	0.08
	Number of sexual partners of susceptible n1 with asymptomatic infective	2.0
	Number of sexual partners of susceptible n2 with symptomatic infective	2.0
	Natural death rate	θ d	0.02
	Removal rate to symptomatic class	α	0.6
	Rate of development to AIDS	κ	0.1

Furthermore, for simplicity, uncertainties of parameters are indicated by uniform distribution in their ranges of variation. The mathematical expectations and uncertainties of input parameters and output R 0 are presented in table 3.6. Regarding the uncertainty in R 0 , both independent and correlated situations are discussed. The underlying sensitivity analysis results are displayed in table 3.7. In our analysis, the first-order Taylor expansions are considered only, which explain 99.6% and 98.8% of the exact uncertainty

Table 3 . 6 :

 36 Uncertainty determination for both input and output parameters. The uncertainty in R0 is obtained by considering the first-order Taylor expansion, which explains 99.6% and 98.8% of the exact uncertainty of R0 for independent and correlated situations, respectively.

		parameter µ (baseline value) Range of variation σ (uncertainty)
		Q0	0.029	[0.0261, 0.0319]	0.002
		β0	0.03	[0.027, 0.033]	0.002
		γ	0.4	[0.36, 0.44]	0.023
		β1	0.2	[0.18, 0.22]	0.012
	input	β2 n1	0.08 2	[0.072, 0.088] [1.8, 2.2]	0.005 0.115
		n2	2	[1.8, 2.2]	0.115
		θ d	0.02	[0.018, 0.022]	0.001
		α	0.6	[0.54, 0.66]	0.035
		κ	0.1	[0.09, 0.11]	0.006
	output	R0	1.429 1.432	--	0.227(ρ = 0) 0.252
					(ρ β 1 n 1

  14, s 2 =0.14, s 3 =0.36, s 12 =0.05, s 13 =0.13, s 23 =0.13, s 123 =0.04} which are in remarkable agreement with Sobol's sampling-

	based values.

The information provided by item [∂y/∂x t ] in Eqs. (4.5) and (4.6) can be approximated

Table 4 . 1 :

 41 Analytic values for uncertainty and sensitivity analysis to numerical test examples. Numbers in parentheses denote the thousandth and ten-thousandth digits of corresponding items.

	y	s 1	s 2	s 3	s T 1	s T 2	s T 3	s 12	s 13	s 23	s 123
	y 1	0.07(14)	0.28(57)	0.64(29)	0.07(14)	0.28(57)	0.64(29)	0	0	0	0
	y 2	0.51(73)	0.31(03)	0.12(93)	0.56(04)	0.31(03)	0.17(24)	0	0.04(31)	0	0
	y 3	0.08(52)	0.08(52)	0.36(37)	0.36(36)	0.36(36)	0.80(11)	0.02(84)	0.18(75)	0.18(75)	0.06(25)
	y 4	0.08(00)	0.08(00)	0.41(57)	0.32(25)	0.32(25)	0.82(15)	0.01(85)	0.18(18)	0.18(18)	0.04(21)

Table 4 . 2 :

 42 Sampling-based results for uncertainty and sensitivity analysis of numerical test examples with M = 10000. Rows of y r 3 and y r 4 indicate the values after introducing our proposed modification. Values provided by LHS-based method are averaged over 1000 simulations with the corresponding standard deviation displayed in parentheses. For the results given by QMC sampling-based method, data in parentheses denote the thousandth and tenthousandth digits.

	LHS-based method

Table 4 . 3 :

 43 Sampling-based values for uncertainty and sensitivity analysis of both SIR and SIS models. The sample size is M = 10000. Rows of SIR r and SIS r indicate the results after introducing our proposed modifications. For LHS-based results, numbers in parentheses denote the standard deviation of each corresponding item, resulting from 1000 independent simulations. For QMC sampling-based values, data in parentheses are the thousandth and ten-thousandth digits. SIR =0.24(.00), E(y) SIS =0.75(.00); V (y) SIR =0.07(.00), V (y) SIS =0.10(.00).

	LHS-based method

SIR =0.24

[START_REF] Borgonovo | Sensitivity analysis: A review of recent advances[END_REF]

, E(y) SIS =0.75(01); V (y) SIR =0.07

[START_REF] Blower | A tale of two futures: Hiv and antiretroviral therapy in san francisco[END_REF]

, V (y) SIS =0.10(42).

Table 5 . 1 :

 51 Sampling-based results for the uncertainty and sensitivity analysis of the model under discussion. The number of samples used here is 10000. Each value offered by LHS-based method is averaged over 100 samplings with number in parentheses denoting the standard deviation resulted from independent samplings. For each sample, the opinion dynamics is established by averaging over 20 random configurations. For QMC method, results are obtained by a deterministic sampling. And the opinion dynamics is established by averaging over 100 random configurations for each sample.

	LHS-based method

results are averaged over 100 samplings used for eliminating the randomness in sample generation. For each sample, the opinion dynamics is established by averaging over 20 random configurations. QMC sampling-based values are provided by a deterministic sampling. And for each sample, the opinion dynamics is established over 100 random configurations . Regarding observation Ω s , its uncertainty is mainly contributed by

  + a 13 σ 2 σ 3 ρ 23 + a 14 σ 2 σ 4 ρ 24 ), (a 12 σ 2 σ 3 ρ 23 + a 13 σ 2 3 + a 14 σ 3 σ 4 ρ 34 ), (a 12 σ 2 σ 4 ρ 24 + a 13 σ 3 σ 4 ρ 34 + a 14 σ 2 -(ρ 13 ρ 23 + ρ 14 ρ 24 ) + (ρ 13 ρ 24 + ρ 14 ρ 23 )ρ 34 -(ρ 12 ρ 23 + ρ 14 ρ 34 ) + (ρ 12 ρ 34 + ρ 14 ρ 23 )ρ 24 (ρ 12 ρ 24 + ρ 13 ρ 34 ) + (ρ 12 ρ 34 + ρ 13 ρ 24 )ρ 23 Similarly, c 11 can be obtained by the definition of the variance of x 1 : + 2a 12 a 13 σ 2 σ 3 ρ 23 + 2a 12 a 14 σ 2 σ 4 ρ 24 + 2a 13 a 14 σ 3 σ 4 ρ 34 ), (B.19) where σ 2 (r 1 ) = σ 2 1 was used. Inserting Eq. (B.17) into Eq. (B.19) provides -ρ 2 23 ) + 2ρ 12 ρ 13 (ρ 23 -ρ 24 ρ 34 ) +2ρ 13 ρ 14 (ρ 34 -ρ 23 ρ 24 ) + 2ρ 12 ρ 14 (ρ 24 -ρ 23 ρ 34 )] 1/2 . (B.20) Through the above analysis of simple cases, general expressions of the coefficients (entries of matrices A and C) that specify the independence and correlation sections separated from arbitrary variables are then derived, see Eqs. (3.16) and (3.17). of x 1 into Eq. (C.3) yields V 1 = 4E [ (a 12 x 2 + a 13 x 3 + c 11 r 1 ) 2 ] = 4(a 2 12 + a 2 13 + 2a 12 a 13 ρ 23 + c 2

	which implies
	c 2 11 = 1 -4 c 11 = 1 σ 2 1 (a 2 12 σ 2 2 + a 2 13 σ 2 3 + a 2 14 σ 2 [ 1 -ρ 2 23 -ρ 2 24 -ρ 2 34 + 2ρ 23 ρ 24 ρ 34 -ρ 2 12 (1 -ρ 2 34 ) -ρ 2 13 (1 -ρ 2 24 ) -ρ 2 ] -1/2 [ = 4(1 -c 2 11 ) + 4c 2 11 , 1 -ρ 2 23 -ρ 2 24 -ρ 2 34 + 2ρ 23 ρ 24 ρ 34 14 (1 11 )	(C.7)
		ρ 12 = 2 ρ 13 = 1 σ 1 σ 2 (a 12 σ 2 1 σ 1 σ 3
		ρ 14 =	1 σ 1 σ 4			4 ),	(B.16)
	which drive
	a 12 = a 13 = a 14 =	ρ 12 (1 -ρ 2 34 ) 1 -ρ 2 23 -ρ 2 24 -ρ 2 34 + 2ρ 23 ρ 24 ρ 34 ρ 13 (1 -ρ 2 24 ) 1 -ρ 2 23 -ρ 2 34 + 2ρ 23 ρ 24 ρ 34 24 -ρ 2 ρ 14 (1 -ρ 2 23 ) -1 -ρ 2 23 -ρ 2 24 -ρ 2 34 + 2ρ 23 ρ 24 ρ 34	σ 1 σ 2 σ 1 σ 3 σ 1 σ 4	, , .	(B.17)
			σ 2 1 = E	[	(x 1 -µ 1 ) 2 ]	,	(B.18)

B.15) Substituting the formulation of x 1 (Eq. (B.14)) into the previous equation, we get

  Analogously, x 2 can be reformed as below based on x 1 and x 3 for the quantification of independent, correlated and coupling variance contributions contained in V 2 :x 2 = a 21 x 1 + a 23 x 3 + c 22 r 2 , (C.9) with r 2 being the element of r that satisfies Eqs. (3.14) and (3.15). Coefficients {a 21 , a 23 , c 22 } determined by Eqs. (3.16) and (3.17). Inserting Eq. (C.9) into Eq. (C.4) drives V 2 =E [ (a 21 x 1 + a 23 x 3 + c 22 r 2 ) 4 ] -E 2 [ (a 21 x 1 + a 23 x 3 + c 22 r 2 ) 2 ] + 2a 21 a 23 x 1 x 3 + 2a 21 c 22 x 1 r 2 + 2a 23 c 22 x 3 r 2 + c 2 22 r 2 2 ) 2 ] -E 2 [ (a 21 x 1 + a 23 x 3 + c 22 r 2 ) 2 ]

	2 11 ,	V C 1 = 4(1 -c 2 11 ),	V UC 1	= 0.	(C.8)
	=E 3 = [ (a 2 21 x 2 1 + a 2 23 x 2 [ 3(a 4 21 + a 4 23 ) + 6a 2 21 a 2 23 (2ρ 2 13 + 1) + 12a 21 a 23 ρ 13 (a 2 21 + a 2 23 )	]
	+ =2(1 -c 2 [ 6c 2 22 (a 2 21 + a 2 23 + 2a 21 a 23 ρ 13 ) 22 ) 2 + 2c 4 22 + 4c 2 22 ), ] + 3c 4 22 -22 (1 -c 2	[	(a 2 21 + a 2 23 + 2a 21 a 23 ρ 13 ) + c 2 22 (C.10) ] 2

  2 (a 12 x 2 + a 13 x 3 + c 11 r1 ) 4 ] + 16E [ x 2 (a 12 x 2 + a 13 x 3 + c 11 r 1 ) 3 ] + 2a 12 a 13 x 2 x 3 + 2a 12 c 11 x 2 r 1 + 2a 13 c 11 x 3 r 1 + c 2 + 2a 12 a 13 x 2 x 3 )r 1 + 3c 2 11 (a 12 x 2 + a 13 x 3 )r 2 1 + c 3 11 r 3 1) the first fraction to the third one, are separately produced by the correlations of x 1 with the rest, the independence of x 1 , and the coupling effect between correlations

	=16E	[	x 2 2 (a 2 12 x 2 2 + a 2 13 x 2 3 11 r 2 1 ) 2 ]
	+ 16E	[	x 2 (a 3 12 x 3 2 + a 3 13 x 3 3 + 3a 2 12 a 13 x 2 2 x 3 + 3a 12 a 2 13 x 2 x 2 3 )	]
	+ 16E	[	x 2 (3c 11 (a 2 12 x 2 2 + a 2 13 x 2 3 ]
	=48(1 -c 2 11 )(1 -c 2 11 + ρ 12 + 4ρ 2 12 ) + 48c 4 11 + 48c 2 11 (2 -2c 2 11 + ρ 12 + 4ρ 2 12 ),	(C.15)
	which, from and independence, specifying
					V C 1 12 = 48(1 -c 2 11 )(1 -c 2 11 + ρ 12 + 4ρ 2 12 ), V U 1 12 = 48c 11 1 4 , V UC 1 12 = 48c 2 11 (2 -2c 2 11 + ρ 12 + 4ρ 2 12 ).	(C.16)

  1 (a 21 x 1 + a 23 x 3 + c 22 r 2 ) 2 ] + 16E [ x31 (a 21 x 1 + a 23 x 3 + c 22 r 2 ) + 30a 21 a 23 ρ 13 + 3a 21 + 3a 23 ρ 13

						]
	= 16	[	15a 2 21 + a 2 23 (12ρ 2 13 + 3) ]	+ 48c 2 22	(C.17)
	= 48(1 -c 2 22 + ρ 12 + 4ρ 2 12 ) + 48c 2 22 ,
	which shows a vanishing coupling variance contribution but the existent correlated (first
	term with brackets) and independent (last term) ones:
	V U 2 12 = 48c 2 22 ,	V C 2 12 = 48(1 -c 2 22 + ρ 12 + 4ρ 2 12 ),	V UC 2 12	= 0.	(C.18)

  3 • (a 12 x 2 + a 13 x 3 + c 11 r 1 ) 2 ] -c 2 11 + 2a 12 ρ 23 (a 13 + a 12 ρ 23 ) + 2a 13 (a 13 + a 12 ρ 23 ) -ρ 2 13 + c 2

	= a 2 12	(	2ρ 2 23 + 1	)	13 -ρ 2 13 + 6a 12 a 13 ρ 23 -ρ 2 + 3a 2 13 + c 2 11	(C.19)
	= 1 11
	= 1 -c 2 11 + ρ 2 13 + c 2 11 ,

  2,1 (x 1 , x 2 2 , x 3 ) -M 2 (x 2 )ρ 13 .By introducing the reforming of x 2 (Eq. (C.9)), we getV 123 =E [ x 1 x 3 (a 21 x 1 + a 23 x 3 + c 22 r 2 ) 2 ] -E [ (a 21 x 1 + a 23 x 3 + c 22 r 2 ) 2 ] ρ 13 =E [ x 1 x 3 (a 21 x 1 + a 23 x 3 ) 2 ] -E [ (a 21 x 1 + a 23 x 3 ) 2 ]

						(C.22)
					ρ 13 + E	[	c 2 22 x 1 x 3 r 2 2	]	(C.23)
	-E	[	c 2 22 r 2 2	]	ρ 13 ,

  a 12 x 2 + a 13 x 3 + a 14 x 4 + c 11 r 1 , (D.2) where r 1 is the element of r that satisfies Eqs.(3.14) and(3.15). Coefficients {a 12 , a 13 , a 14 , c 11 } are presented in Eqs. (B.17) and (B.20). Inserting (D.2) into (D.1) providesM 1,1,1,1 (x 1 , x 2 , x 3 , x 4 ) = E [ (a 12 x ′ 2 + a 13 x ′ 3 + a 14 x ′ 4 + c 11 r ′ = x i -µ i for i = 2, 3, 4 and r ′ 1 = r 1 -µ(r 1 ). In the determination of E [ a 12 (x ′ 2 ) 2 x ′ 3 x ′ 4 ], another variable should also be reformed based on the others enclosed in the brackets. Here we formulate x 3 (equivalent to formulate x 2 or x 4 as input variables are normally distributed and hold the same standard deviation) asx 3 =a 32 x 2 + a 34 x 4 + c 33 r 3 , (D.4) where r 3 is the element of r that satisfies Eqs. (3.14) and (3.15). Coefficients {a 32 , a 34 , c 33 } are determined via Eqs. (3.16) and (3.17) as a 32 = (ρ 23 -ρ 24 ρ 34 )/(1 -ρ 2 24 ), (D.5) a 34 = (ρ 34 -ρ 23 ρ 24 )/(1 -ρ 2 24 ), (D.6) c 33 = (1 -ρ 2 23 -ρ 2 24 -ρ 2 34 + 2ρ 23 ρ 24 ρ 34 ) 1/2 /(1 -ρ 2 24 ) 1/2 . 32 (x ′ 2 ) 3 x ′ 4 + a 12 a 34 (x ′ 2 ) 2 (x ′ 4 ) 2 ] = 2a 12 ρ 24 (a 32 + a 34 ρ 24 ) + a 12 (a 34 + a 32 ρ 24 ) = 2a 12 ρ 23 ρ 24 + a 12 ρ 34 .The other two average items in Eq. (D.3) can be similarly determined as ] = 2a 13 ρ 23 ρ 34 + a 13 ρ 24 , (D.9) ] = 2a 14 ρ 24 ρ 34 + a 14 ρ 23 .(D.10)We now getM 1,1,1,1 (x 1 , x 2 , x3 , x 4 ) =2(a 12 ρ 23 ρ 24 + a 13 ρ 23 ρ 34 + a 14 ρ 24 ρ 34 ) + (a 12 ρ 34 + a 13 ρ 24 + a 14 ρ 23 ) =ρ 34 (a 12 + a 13 ρ 23 + a 14 ρ 24 ) + ρ 24 (a 13 + a 12 ρ 23 + a 14 ρ 34 ) + ρ 23 (a 14 + a 12 ρ 24 + a 13 ρ 34 ) =ρ 12 ρ 34 + ρ 13 ρ 24 + ρ 14 ρ 23 .Analogously, the main variance contribution produced by x 2 is represented asV 2 = • [M i+j (x 2 ) -M i (x 2 )M j (x 2 )] , (E.5)which is provided by the second term of Eq. (E.1): 7 sin 2 (x 2 ). Zero mean value of x 2

	∞ ∑ i,j=0 )({µ}) suggests 1 i! • j! ( ∂ i y ∂x i 2 • ∂ j y ∂x j 2
	(	∂ i sin 2 (x 2 )
		∂x i 2	µ 2 =0
	1 )x ′ 2 x ′ 3 x ′ 4 3 ) 2 x ′ 2 (x ′ a 13 x ′ 4 ] + E a 12 a (D.8) ] a 12 (x ′ 2 ) 2 x ′ 3 x ′ 4 ] + E [ [ a 14 x ′ 2 x ′ 3 (x ′ 4 ) 2 ] , (D.3) = [ a 13 x ′ = E [ 2 ) 2 x ′ a 12 (x ′ Then we can obtain where x ′ E [ 3 x ′ 4 ] E [ 2 (x ′ 3 ) 2 x ′ 4 i (D.7) E [ a 14 x ′ 2 x ′ 3 (x ′ 4 ) 2 (D.11)

  1 • π 4 (1 -ρ 2 13 ) 4

				9(2i + 2j + 3)	+	2(1 -ρ 2 13 ) 2 5(2i + 2j + 3)	]	,	(E.26)
	V C 3 13 = 0.1	∞ ∑ i,j=0	(-1) i+j π 2(i+j+1)+4 (2i + 1)!(2j + 1)!	[ 0.1 • π 4 ρ 8 13 2i + 2j + 11	+	2ρ 4 13 2i + 2j + 7	]	.	(E.27)

  F.2. The fraction of agents who will migrate in opinion from R2 and R3 to R1 is the same as that who will migrate from R2 to R1 in the previous

	situation, that is	n 2,3→1 =	1 2	ε 2 .	(F.10)
	This implies the fraction of winners at steady state is given by
		n w =	3 8	+	1 2	ε 2 .	(F.11)

  Similarly, denote the maximal opinion among the group of winners as s. The final average opinion in clique G -then can be represented by

										4	) 2 ].	(F.14)
	Now we get the final average opinion in clique G + :
	Ω + s = Ω 3 s + Ω 4 s = -	1 6	(ε -	1 4	) 3 +	1 4	.	(F.15)
		Ω -s = n w •	1 2	(s -1) + Ω 2 s .	(F.16)
	By recalling Eq. (F.7), we obtain
		Ω -s = 1 -	1 8ε 2 -26	(ε -	391 12	),	(F.17)
	which, together with previously calculated Ω + s , specifies
	Ω s =	1 8ε 2 -26	(ε -	391 12	) +	1 6	(ε -	1 4	) 3 -	5 4	.	(F.18)

If the confidence threshold is limited to the interval (1/2, 3/4], we can divide the opinion

  Solution of the above inequality (|s| ≥ 1) implies the clustering of all agents in opinion -1 at steady state. So we get Ω s = 1. (F.31) Résumé Par un développement en série de taylor, une relation analytique générale est établie pour calculer l'incertitude de la réponse du modèle, en assumant l'indépendance des entrées. En utilisant des relations de puissances et exponentielles, il est démontré que l'approximation souvent utilisée permet d'évaluer de manière satisfaisante l'incertitude sur la réponse du modèle pourvu que l'incertitude d'entrée soit négligeable ou que le modèle soit presque linéaire. La méthode est appliquée à l'étude d'un réseau de distribution électrique et à un modèle d'ordre économique. La méthode est étendue aux cas où les variables d'entrée sont corrélées. Avec la méthode généralisée, on peux déterminer si les corrélations d'entrée doivent ou non être considérées pour des applications pratiques. Des exemples numériques montrent l'efficacité et la validation de notre méthode dans l'analyse des modèles tant généraux que spécifiques tels que le modèle déterministe du vih. La méthode est ensuite comparée à celle de sobol. Les résultats montrent que la méthode de sobol peut surévaluer l'incidence des divers facteurs, mais sousestimer ceux de leurs interactions dans le cas d'interactions non linéaires entre les paramètres d'entrée. Une modification est alors introduite, aidant à comprendre la différence entre notre méthode et celle de sobol. Enfin, un modèle numérique est établi dans le cas d'un jeu virtuel prenant en compte la formation de la dynamique de l'opinion publique. L'analyse théorique à l'aide de la méthode de modification d'un paramètre à la fois. La méthode basée sur l'échantillonnage fournit une analyse globale de l'incertitude et de la sensibilité des observations. Mesure de sensibilité, Dynamique d'opinion

		This requires
						|s| 2 (|s| + 1) N • 1	• R -≥	R N	.	(F.30)
	Mots clés Analyse	d'incertitude,	Analyse	de	sensibilité,
	Décomposition	de	variance,	Echantillonnage,
	Corrélation,

Acknowledgements

My Ph.D. study of three years is near its ending. I was continuously supported and encouraged by my supervisors, colleagues, friends, and family members. At this special and important moment, I would like to express my sincere gratitude to all of them.

, I am grateful to Dr. Jian Jiang for his valuable comments in the writing of the thesis. I gratefully acknowledge the financial support that enabled me to finish smoothly my Ph.D. study in France. I was funded by China Scholarship Council. My work was also partially supported by the Programme of Introducing Talents of Discipline to Universities of China, the National Natural Science Foundation of China, the Ecole Doctorale 3M at Le Mans University, and the Institut des Molécules et Matériaux du Mans (IMMM).

Abbreviations vi

Abbreviations AIDS where (x 1 , x 2 , x 3 ) ∼ N (µ, Σ) with mean vector µ = (0, 0, 0) T The variance V (y) of model response can be similarly computed by the use of Eq. (3.4) as V (y) = 55 + 48ρ 12 + 2ρ 12 ρ 23 + 192ρ 2 12 + ρ 2 13 ,

which is generated by the partial variance contributions involving

(3.50)

The independent, correlated and coupling variance contributions divided from the main variance contributions V 1 and V 2 are stated as

).

(3.51)

Regarding the existent higher-order partial variance contributions, we have = 0,

where i ∈ {1, 3}, j ∈ {1, 2, 3} and c i are determined with Eq. (3.17). A detailed calculation process for the above items are presented in Appendix C. In table 3.2, the exact sensitivity indices are listed, showing a dominated influence of the interaction effect between x 1 and x 2 in the absence of input correlations. In the presence of input correlations, the independent, correlated and coupling effects provided by x 1 and x 2 are all significant in establishing the uncertainty of model response. The convergence of our analysis results with the order of Taylor expansion is also presented in Fig. 3.1. γ is the highest-order of Taylor expansion that is considered in our method. Analysis results are convergent at γ = 3 (the highest-order of Taylor expansion of the model output), implying the effects of any possible orders of Taylor expansion of the model output are and

respectively. Simulation results suggest N + /N -= R + /R -at steady state (Fig. 5.3(a)).

Inserting R + = γR/(1 + γ) and R -= R/(1 + γ) into Eqs. (5.15) and (5.16) yields

(5.17)

which determines the order parameter as

(5.18)

Sufficiently biased resource allocation, that is γ < 1/3 obtained from the inequation

will result in the absence of winners in G + at the initial time. However, some agents will be attracted from G + to G -by communicating with their friends. The decrease in the number of agents with positive opinions will then generate winners in G + . The distribution of the fraction of winners presented in Fig. 5.3(b) demonstrates that almost all agents are winner at final state when 10/N < γ < 1/3 (γ = 0 is also included).

Assume agents in G + gather finally in opinion s. The condition for winning requires This implies the final opinions in R3 is non-uniformly distributed and

(5.36)

The sum of Ω 3 s and Ω 4 s yields the final average opinion in G + ,

(5.37)

Then we have

(5.38)

For any other possible values of ε, Ω s can be analogously determined and is analysed in Appendix F in detail. Finally, we obtain

(5.39) Similar analysis is also performed for scale-free networks with k a = 32 (as an example of dense networks) and fully connected networks, that is, k a = N -1. When k a = 32, the order parameter is expressed as

(5.40) present, nevertheless, sampling-based strategies regarding uncertainty and sensitivity analysis have been designed by assuming that input variables are independent of each other. Following the presence of input correlations in practical applications, we shall aim at extending the sampling-based strategy to the case of input correlations. This may contribute to the application of uncertainty and sensitivity analysis in the modeling specific to more general practical problems, and help decide whether or not the input correlations should be considered in practice when the functional form connecting input and output variables is absent.

Sensitivity analysis for the economic time series

For now, the uncertainty and sensitivity analysis is mainly implemented for models with scalar output. In many situations, however, the output variable is not always a scalar but could also be a function. This commonly happens in time-dependent systems. A most typical example of time-dependent systems is the economic time series. In official statistics and applied macroeconomics, it has been customary to analyse an economic time series by extracting from it a long term movement, or trend, for separate study and then scrutinizing the residual portion for short term oscillatory movements and random fluctuations. With given time series, we shall try to quantify the roles of some possible parameters in the generation of a typical shape and fluctuations, and in trend extraction, using Sobol's measures.

Sensitivity analysis for the spread of an epidemic disease

The spread of an epidemic disease is referred to as another typical time-dependent problem. In the context of uncertainty and sensitivity analysis of models for the reproduction of disease spreading, what researchers are usually interested in is the equilibrium state of the underlying system, as that we discussed in Chapter 4. Focusing on the equilibrium state of epidemic spreading neglects the influence of temporal dimension. In practical applications, however, the evolution process of epidemic disease with time parameter is much more important than the equilibrium state, which provides fruitful information for identifying the outbreak of a disease and then makes its efficient controlling possible.

We shall attempt to quantitatively identify the global roles of different parameters in the propagation process of an epidemic disease along the axis of time, whereby the outbreak of a disease could be understood more deeply.

Appendix A

Central moments

The kth central moment of a continuously distributed variable x is defined as

(A.1)

A.1 Uniform distribution

Suppose x is uniformly distributed. Its PDF is then given by

which provides µ = (x m + x 0 )/2, and standard deviation σ = (x m -x 0 )/(2 √ 3). Substituting Eq. (A.2) into Eq. (A.1) yields

] ,

which can be simplified as follows by inserting µ = (x m + x 0 )/2:

The above expression suggests M k (x) keeps at zero for odd k but acts as a function of σ for even k, that is

k+1 σ k , k is even.

(A.5)

A.2 Normal distribution

Regarding normally distributed x, the PDF is defined as

By inserting the above expression into Eq. (A.1), we have

which fixes at zero for odd k since it involves the integral of an odd function over a real line. While k is even, M k (x) is obtained as

where k!! denotes the double factorial of k with 0!! = 1. The kth central moment of a normally distributed variable x is obtained as

Appendix B

Generation of correlated variables B.1 Two correlated variables

Suppose x 1 and x 2 are only correlated with each other in the input space with given correlation coefficient ρ 12 . By using the linear correlation model, x 1 can be formulated on the basis of x 2 as

where r 1 is the element of r that satisfies Eqs. (3.14) and (3.15). Employing the definition of correlation coefficient, see Eq. (3.9), we have

which yields

Furthermore, we also have

By substituting σ 2 1 = V (r 1 ), coefficient c 11 is obtained as

Derivation of Eqs. (3.51-3.52)

The form of related model function is

where (x 1 , x 2 , x 3 ) ∼ N (µ, Σ) with mean vector µ = (0, 0, 0) T and covariance matrix

C.1 First-order contributions

With Eq. (3.4), fractional variance contribution produced by each input factor alone are represented as

2

To determine the independent, correlated and coupling variance contributions that are contained in V 1 , x 1 should be reformed on the bases of x 2 and x 3 as

in which r 1 is the element of r that satisfies Eqs.

Ishigami function

The Ishigami function is expressed as

with input variables uniformly distributed in the interval [-π, π] which provides µ i =0 and σ i = π/ √ 3 for i = 1, 2, 3.

E.1 First-order contributions

With help of Eq. (3.4), the main partial variance contribution produced by input x 1 alone is represented as

The existent partial derivatives of y with respect to x 1 are totally provided by the first term of Eq. (E.1): sin(x 1 ) as µ 3 = 0. Zero mean value of x 1 determines

Substituting Eqs. (A.5) and (E.3) into Eq. (E.2) yields

The independent and correlated variance contributions are then evaluated by setting k = 0 and k = i + j + 1, respectively:

The coupling variance contribution is then evaluated accordingly by

(E.13)

E.2 Second-order contributions

For partial variance contributions of second-order, the form of Ishigami function provides

3

3

3 )Mj,4(x j 1 , x 4

3 ) ]

3 ) -Mi(x1)Mj,4(x j 1 , x 4

3 ) ]

By substituting Eq. (E.3), V 13 is updated as

To determine V 13 , it is necessary to derive the covariance items enclosed in the above square parentheses.

which determines the final fraction of agents in R1 as

Analogous to Eq. (5.36), at steady state, the average opinion in the region R2 can be specified as

With respect to the limitation of confidence threshold, opinion distributions in both R3 and R4 will keep at their initial configurations (uniformly distributed) throughout the space into five adjacent regions (Fig. F.3). As stated before, all agents holding initial opinions in R2 will migrate to R1 in opinion by communicating with their friends. And the fraction of agents who will be affected in opinion by their friends from R3 and R4 to R1 is approximated as

)dx = 1 8 .

(F.19)

The fraction of agents holding steady opinions in R1 is then determined as

The steady average opinion in G + is the same as that in the previous case, i.e., Eq.

(F.15), since the partition of opinion space in the positive interval keeps unchanged. For region R3, its average opinion at steady state is analogously computed as

(F.21)

Agents taking opinions in R1 will finally win at gambling by further local communications. This suggests the average opinion in G -will stabilize at For high confidence threshold: ε ∈ ( 5 4 , 7 4 ], the opinion space can be divided into three regions (Fig. F.4). The fraction of agents who will be affected in opinion by their friends