La technique développée par le mathématicien Russe I. M. Sobol va permettre de prendre en compte toutes les interactions entre les n grandeurs d’entrée, en n’ayant aucune connaissance sur l’expression du modèle de mesure f. Ainsi nous allons obtenir l’indice de sensibilité total pour chaque grandeur d’entrée. Sobol s’est inspiré de divers travaux de statistique [1, 2] et les a généralisés pour fournir une mise en oeuvre directe du concept de la méthode Monte Carlo [3].

Pour récapituler, un mesurande Y avec pour modèle mathématique f est donnée par :

avec X_i(i = 1, ..., n) ses paramètres d’entrée. Sobol considère la décomposition fonctionnelle de f pouvant s’exprimer comme une somme de fonctions de dimension croissante [3], où :

avec :

L’espérance 𝔼(Y|X_i = q_i) est la moyenne de Y où seulement les valeurs de q_i engendrées par l’espace de la i-ème grandeur d’entrée X_i considérée sont conditionnées (fixées). De même que pour l’espérance 𝔼(Y|X_i,X_j) où seulement les valeurs de q_i et de q_j des grandeurs d’entrée X_i et X_j considérées sont conditionnées (fixées), et ainsi de suite. Cette décomposition est appelée "High Dimensional Model Representations (HDMR) [4]. Sobol impose certaines contraintes [3] afin d’avoir l’unicité de cette décomposition et par conséquent d’assurer la décomposition de la variance globale en une somme de variances partielles :

avec :

Avec V[𝔼(Y|X_i)] qui est la variance de Y en ayant maintenu fixe les valeurs de X_i considérée (i-ème grandeur d’entrée). De même pour V[𝔼(Y|X_i,X_j)] qui est la variance de Y en ayant maintenu fixe les valeurs de X_i et X_j considérées. Ensuite, en divisant tous les termes par la variance du mesurande V(Y) :

Par conséquent, Sobol se base sur cette décomposition pour définir tous les indices de sensibilité d’ordre 1 et supérieur de toutes les grandeurs d’entrée X_i :

Ainsi, l’indice de Sobol de premier ordre mesure l’effet principal pour une grandeur d’entrée X_i considérée :

Puis toutes ses interactions avec les autres grandeurs d’entrée, telles que l’ordre 2 entre X_i et X_j considérées :

La subtilité pour calculer la variance V[𝔼(Y|X_i)] est de construire, par la méthode de Monte Carlo, deux matrices indépendantes [^①q₁^r, ^①q₂^r, ..., ^①q_n^r] et [^❷q₁^r, ^❷q₂^r, ..., ^❷q_n^r] qui constituent deux échantillonages des n grandeurs d’entrée X_i [4, 5, 6]. Ainsi, pour ne prendre en compte que les variations sur les valeurs q_i de la i-ème grandeur X_i considérée, seulement les valeurs de cette grandeur d’entrée seront conditionnées, c’est à dire que nous fixons les valeurs de X_i (grandeur d’entrée considérée) grâce à la première matrice mais que les valeurs des autres grandeurs d’entrée vont être rééchantillonnées en utilisant la deuxième matrice. Pour cela, nous construisons les mesurandes ^①Y , ^❷Y et ^①❷Y en conditionnant une grandeur d’entrée X_i considérée :

Ainsi, la méthode de Sobol est de réécrire l’indice de sensibilité de premier ordre Si de la manière suivante :

où la variance se calculera grâce aux mesurandes ^①Y et ^①❷Y. En remplaçant la variances et covariance dans 11 avec les nouveaux mesurandes ^①Y et ^①❷Y, nous obtenons un estimateur empirique de Sobol de l’indice de sensibilité S_i :

Les moyennes des deux mesurandes ^①Y^r et ^①❷Y_i^r sur la totalité des run r = 1, ..., M sont pratiquement les mêmes, ainsi 𝔼(^①Y) = 𝔼(^①❷Y). Finalement, l’estimateur empirique proposé par Sobol pour calculer les indices de sensibilité de premier ordre est donné par :

Ainsi, les valeurs du mesurande ^①Y^r sont multipliées aux valeurs du mesurande ^①❷Y_i^r pour lesquelles toutes les grandeurs d’entrée ont été rééchantillonnés sauf pour les valeurs de la grandeur d’entrée X_i considérée qui ont été prises dans ^①Y^r. Par conséquent, si la grandeur X_i considérée n’a pas d’influence, alors les valeurs élevée de ^①Y^r sont associées de manière aléatoire avec les valeurs hautes de ^①❷Y_i^r, c’est à dire que ce sont les autres grandeurs d’entrée qui imposent les valeurs hautes (positives et négatives) de ^①Y^r et ^①❷Y_i^r. Par contre, si la grandeur X_i considérée a un impact de premier ordre dans le modèle f, les valeurs élevées de ^①Y^r seront multipliées préférentiellement par des valeur élevées de ^①❷Y_i^r, en augmentant ainsi la valeur de leur produit résultant (^①Y^r · ^①❷Y_i^r). C’est à dire que cette fois-ci c’est X_i considérée qui impose les valeurs hautes (positives et négatives) de ^①Y^r et ^①❷Y_i^r.

Cette méthode permet donc de mesurer les indices de Sobol de premier ordre pour chaque grandeur d’entrée X_i en conditionnant à chaque fois la grandeur d’entrée X_i choisie. Pour mesurer les indices supérieurs de sensibilité, il suffit de la même façon de conditionner deux ou plusieurs grandeurs d’entrée pour obtenir leur sensibilité. Cependant, le nombre d’indices de sensibilité est égal à 2ⁿ−1 en prenant en compte toutes les interactions d’ordres 1 à n. Lorsque le nombre de grandeur d’entrée n est important tels que dans nos travaux, le nombre d’indices de sensibilité devient beaucoup trop grand.

Néanmoins, pour limiter le nombre de calcul et obtenir les indices de sensibilité des interactions, une autre méthode consiste à mesurer directement l’indice de sensibilité total S_Ti de la grandeur d’entrée X_i et de la soustraire à l’indice de premier ordre Ŝ_i calculé précédemment [6]. L’indice de sensibilité total comprend à la fois l’indice de sensibilité de premier ordre de X_i (grandeur d’entrée considérée) et toutes ses interactions avec les autres grandeurs d’entrée.

A partir de l’identité connue [7] où la variance est décomposée entre l’effet principal et résiduel, nous avons :

Le terme V(𝔼(Y|X_i)) mesure l’effet d’ordre 1 de la grandeur d’entrée X_i, alors que le terme 𝔼(V(Y|X_i)) est habituellement appelé le résidu. Cette fois-ci, au lieu de ne conditionner qu’une seule grandeur d’entrée X_i, nous allons conditionner toutes les grandeurs d’entrée, sauf une X_~i [8, 9] :

Le terme V(𝔼(Y|X_~i)) est l’effet d’ordre un de X_~i, c’est à dire que ce terme est la variance de l’espérance de Y conditionnant toutes les variables sauf celles de X_i considérée. Ainsi V(Y) − V(𝔼(Y|X_~i)) enlève la contribution de tous les termes dans la décomposition de variance (equation 4) qui incluent toutes les grandeurs d’entrée du modèle sauf les termes où intervient la grandeur X_i considérée. Ainsi, le terme résiduel 𝔼(V(Y|X_~i)) devient donc la variance restante de tous les termes où intervient X_i considéré (ordre 1 et interactions). En divisant par V(Y) nous obtenons l’indice total de sensibilité :

S_Ti mesure l’effet total qui comprend les sensibilités d’ordre 1 et supérieur (interactions) de la grandeur X_i considérée. La méthode est la même pour calculer V(𝔼(Y|X_~i)) en utilisant les deux matrices ^❷Y^r et ^①❷Y_i^r construites précédemment, où la variance est conditionnée cette fois-ci par rapport à tous les grandeurs d’entrée, sauf une X_~i, c’est à dire que cette fois-ci, seulement la grandeur d’entrée X_i considérée sera rééchantillonée.

De la même façon que précédemment, l’indice total de sensibilité S_Ti peut se réécrire de la manière suivante :

En remplaçant la variance et covariance dans 17 avec les nouveaux mesurandes ^❷Y et ^①❷Y, nous obtenons un estimateur empirique de l’indice de sensibilité totale S_Ti :

Ainsi, les valeurs du mesurande ^❷Y^r sont multipliées aux valeurs du mesurande ^①❷Y_i^r pour lesquelles toutes les grandeurs d’entrée ont été figées (prises dans ^❷Y^r) sauf pour les valeurs de X (grandeur d’entrée considérée) qui ont été rééchantillonnées (prises dans ^①Y^r).

Par définition S_Ti est plus grand que Si puisqu’il intègre à la fois l’effet principal et toutes les interactions de la grandeur d’entrée X_i. Par conséquent, si nous souhaitons seulement l’indice de sensibilité des interactions, il suffit de soustraire S_Ti avec S_i. Par ailleurs, la somme des indices totaux de chaque grandeur d’entrée ne fait pas 1 car par exemple la sensibilité de l’interaction entre trois paramètres S_1,2,3(X₁,X₂,X₃) va être prise en compte à la fois dans l’indice total S_T1 , S_T2 et S_T3 . Dans le cas où la grandeur d’entrée n’a aucune influence sur le modèle, son indice total de sensibilité est nul, elle peut donc être fixée à une constante durant la construction des deux matrices, sans affecter la variance du mesurande pour le calcul de sensibilité des autres grandeurs d’entrée. C’est pourquoi nous utilisons au préalable le plan de Morris afin de fixer constantes toutes les grandeurs d’entrée non influentes. L’intérêt est de réduire le nombre de calcul pour calculer l’effet total de chaque X_i. En effet, par exemple dans le cas où nous avons n = 100 grandeurs d’entrée et nous souhaitons connaitre leur indice de sensibilité total. Nous construisons au préalable deux matrices avec chacune M = 1000 run :

A partir de ces deux matrices, les mesurandes ^①Y^r et ^❷Y^r peuvent ainsi être construits (avec r = 1, ..., M). Il reste maintenant à mettre en place tous les mesurandes ^①❷Y_i^r (avec i = 1, ..., 100 et r = 1, ..., M) afin de pouvoir calculer l’indice de sensibilité total des n = 100 grandeurs d’entrée :

Nous avons donc en tout 102000 run à calculer pour déterminer les valeurs des mesurandes ^①❷Y_i^r. Une fois les mesurandes en place, il suffit d’appliquer l’équation 18 pour obtenir l’indice total de sensibilité de chaque grandeur d’entrée. Dans le cas où nous souhaitons une meilleure estimation de l’indice de sensibilité, cette opération peut être répétée plusieurs fois (permet de converger l’indice de sensibilité et de l’associer à un écart-type) et le nombre de run à calculer peut très vite monter à 10⁶ dans le cas de 10 répétitions. D’où la nécessité d’utiliser au préalable le plan de Morris pour éviter une surcharge de calcul.

En conclusion, le principe de l’estimation des indices de sensibilité, selon la méthode de Sobol, est de simuler deux échantillons à partir des grandeurs d’entrée, puis d’estimer la variance conditionnelle par rapport à chaque grandeur d’entrée considérée en combinant ces deux échantillons. Le principal inconvénient de la méthode de Sobol est son temps de calcul puisqu’elle nécessite une taille d’échantillon suffisamment grande pour assurer la convergence des indices selon la complexité du modèle de mesure. Cependant, cet inconvénient est rapidement contourné par l’utilisation du plan de Morris. Finnallement, le grand avantage de l’analyse de sensibilité avec la méthode de Sobol est qu’elle permet de ne faire aucune hypothèse sur les propriétés de la fonction f du modèle mathématique, elle est adaptée quelque soit le modèle (non-linéaire, non-monotone, code de calcul, ...).

[1] : Sacks J, Welch W J, Mitchell T J and Wynn H P 1989 Design and analysis of computer experiments Statistical Science 4 409–423
[2] : Cukier R, Levine H and Shuler K 1978 Nonlinear sensitivity analysis of multiparameter model systems Journal of Computational Physics 26 1–42
[3] : Sobol' I 1993 Sensitivity estimates for nonlinear mathematical models Mathematical Modelling and Computational Experiments 1 407–414
[4] : Sobol' I 2001 Global sensitivity indices for nonlinear mathematical models and their monte carlo estimates Mathematics and Computers in Simulation 55 271–280
[5] : Saltelli A 2002 Making best use of model evaluations to compute sensitivity indices Computer Physics Communications 145 280–297
[6] : Homma T and Saltelli A 1996 Importance measures in global sensitivity analysis of nonlinear models Reliability Engineering & System Safety 52 1–17
[7] :Mood A M, Graybill F A and Boes D C 1974 Introduction to the Theory of Statistics 3rd Edition (Tmhe Ie Overruns) (McGraw Hill) ISBN 0070854653
[8] : Jansen M J 1999 Analysis of variance designs for model output Computer Physics Communications 117 35–43
[9] : Saltelli A, Annoni P, Azzini I, Campolongo F, Ratto M and Tarantola S 2010 Variance based sensitivity analysis of model output. design and estimator for the total sensitivity index Computer Physics Communications 181 259–270