[PYTHON] Analyse de sensibilité globale à l'aide des bibliothèques d'analyse de sensibilité Salib et Oacis

Qu'est-ce que l'analyse de sensibilité?

L'analyse de sensibilité est une méthode pour évaluer quantitativement quelle incertitude d'entrée est dérivée de l'incertitude de la sortie d'un modèle ou d'un système mathématique. https://en.wikipedia.org/wiki/Sensitivity_analysis

Lorsque la sortie est $ Y $ et l'entrée est $ \ bf X $, le modèle mathématique est exprimé comme suit. S'il y a de l'incertitude dans une certaine entrée $ X_i $, l'incertitude apparaît également dans $ Y $, mais quantitativement combien l'incertitude de $ Y $ dépend de laquelle $ X_i $ est incertaine. évaluer.

Y = f({\bf X})

Il est utilisé aux fins suivantes.

Méthode d'analyse de sensibilité

Il existe plusieurs méthodes d'analyse de sensibilité. Utilisez correctement en fonction du coût de calcul, de l'objectif et des caractéristiques du modèle mathématique.

Une à la fois (OAT) est la méthode la plus simple. Comme son nom l'indique, il s'agit d'une méthode pour vérifier à quel point le résultat change en ne déplaçant qu'un seul paramètre par une méthode telle que le coefficient différentiel partiel ou la régression linéaire. C'est très simple et facile à comprendre, et le coût de calcul est faible, mais comme les paramètres autres que les paramètres mobiles sont fixés aux valeurs de base, l'espace d'entrée ne peut pas être entièrement recherché.

De plus, lorsque l'on essaie d'évaluer par le coefficient différentiel partiel, la valeur à la ligne de base spécifique ne peut être évaluée que localement. Impossible d'évaluer l'interaction entre les variables d'entrée ou lorsque $ y $ change de manière non monotone par rapport à $ x_i $

global sensitivity analysis (Sobol method)

L'analyse de sensibilité basée sur la variance est une méthode permettant d'effectuer une analyse de sensibilité globale pour résoudre ces problèmes. Elle est souvent appelée la méthode Sobol. https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis

Dans cette méthode, les entrées et les sorties sont considérées comme des variables stochastiques. Décomposez la distribution de la sortie en contributions d'une entrée (ou de plusieurs entrées). Par exemple, on peut dire que 70% de la dispersion de sortie provient de la fluctuation de $ X_1 $, 20% provient de $ X_2 $, et les 10% restants proviennent de l'interaction entre $ X_1 $ et $ X_2 $.

Cette méthode peut gérer (i) des réactions globales et (ii) non linéaires, et (iii) mesurer l'interaction des paramètres d'entrée (effets qui apparaissent lorsque plusieurs paramètres sont modifiés en même temps). Il y a. D'autre part, comme il est traité statistiquement, il est nécessaire d'échantillonner suffisamment l'espace d'entrée et, en général, un grand nombre de calculs avec diverses entrées sont nécessaires et le coût de calcul est lourd.

Quel type de montant doit être calculé spécifiquement? Soit maintenant le modèle une variable de dimension $ d $ $ \ bf X = \ {X_1, X_2, \ dots, X_d } $. On suppose que toutes les variables sont normalisées comme $ X_i \ in [0,1] $ et sont uniformément réparties. La sortie $ Y $ peut généralement être écrite comme suit:

Y=f_{0}+\sum _{i=1}^{d}f_{i}(X_{i})+\sum _{i<j}^{d}f_{ij}(X_{i},X_{j})+\cdots +f_{1,2,\dots ,d}(X_{1},X_{2},\dots ,X_{d})

Où $ f_i $ est une fonction de $ X_i $ et $ f_ {ij} $ est une fonction de $ X_i $, $ X_j $.

Si vous calculez la distribution de Y $

\operatorname {Var}(Y)=\sum _{{i=1}}^{d}V_{i}+\sum _{{i<j}}^{{d}}V_{{ij}}+\cdots +V_{{12\dots d}}

Peut être écrit comme la somme des termes dont la variance dépend de chaque $ i $. Ici, $ V_i $ est défini par la formule suivante.

V_{{i}}=\operatorname {Var}_{{X_{i}}}\left(E_{{{\textbf  {X}}_{{\sim i}}}}(Y\mid X_{{i}})\right),

Cette formule fixe d'abord $ X_i $ à une valeur spécifique et l'intègre à des variables autres que $ i $ pour calculer la valeur attendue de $ Y $ pour un $ X_i $ spécifique, et ce montant de $ X_i $ Cela signifie calculer la variance lorsque le est modifié.

Le terme quadratique $ V_ {ij} $ est également défini comme suit.

V_{ij} = \operatorname{Var}_{X_{ij}} \left( E_{\textbf{X}_{\sim ij}} \left( Y \mid X_i, X_j\right)\right) - \operatorname{V}_{i} - \operatorname{V}_{j}

Il existe principalement deux types d'indices pour évaluer la sensibilité du paramètre $ i $: l'indice du premier ordre $ S_i $ et l'indice d'effet total $ S_ {Ti} $. Chacun est défini comme suit.

S_{i}={\frac  {V_{i}}{\operatorname {Var}(Y)}}

S_{{Ti}}={\frac  {E_{{{\textbf  {X}}_{{\sim i}}}}\left(\operatorname {Var}_{{X_{i}}}(Y\mid {\mathbf  {X}}_{{\sim i}})\right)}{\operatorname {Var}(Y)}}=1-{\frac  {\operatorname {Var}_{{{\textbf  {X}}_{{\sim i}}}}\left(E_{{X_{i}}}(Y\mid {\mathbf  {X}}_{{\sim i}})\right)}{\operatorname {Var}(Y)}}

$ S_i $ est une méthode pour évaluer la distribution de $ Y $ lorsque $ i $ "seulement" est changé. Lorsque $ i $ est modifié seul, il fluctue d'environ $ S_i $ en moyenne. D'autre part, l'indice d'effet total est la contribution de $ i $ lorsque tous les paramètres sont modifiés, y compris la contribution de l'interaction avec d'autres paramètres. Si $ S_ {Ti} $ est petit, la valeur n'affectera pas le résultat même si elle est fixe, elle peut donc être utilisée pour simplifier le modèle.

Pour calculer ces indicateurs, l'espace des paramètres est échantillonné à Monte Carlo. Un échantillonnage aléatoire peut être utilisé, mais la méthode Quasi-Monte Carlo est souvent utilisée comme un dispositif pour un calcul plus efficace. En fait, les bibliothèques décrites ci-dessous le font. Habituellement, un échantillonnage de plusieurs milliers de x (nombre de paramètres) est nécessaire.

Comment utiliser SALib

Une des bibliothèques qui effectue une analyse de sensibilité globale est SALib. Ici, comment l'utiliser sera expliqué. Plusieurs méthodes sont implémentées dans cette bibliothèque, dont la méthode Sobol, mais seule la méthode Sobol sera expliquée ici.

Référence: https://www.youtube.com/watch?v=gkR_lz5OptU&feature=youtu.be

La méthode d'installation est pip install SALib Nécessite numpy, scipy, matplotlib comme prérequis.

SALib est effectué dans les 4 étapes suivantes.

  1. Définition de l'entrée du modèle
  2. Génération d'échantillons
  3. Exécution du modèle mathématique pour lequel une analyse de sensibilité est souhaitée
  4. Calculez les indicateurs d'analyse de sensibilité ($ S_i $ et $ S_ {Ti} $)

L'exemple de code le plus simple est le suivant. Ici, une fonction de test appelée ʻIshigami` est utilisée comme modèle mathématique.

from SALib.sample import saltelli
from SALib.analyze import sobol
from SALib.test_functions import Ishigami
import numpy as np

# Define the model inputs
problem = {
    'num_vars': 3,
    'names': ['x1', 'x2', 'x3'],
    'bounds': [[-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359]]
}

# Generate samples
param_values = saltelli.sample(problem, 1000)

# Run model (example)
Y = Ishigami.evaluate(param_values)

# Perform analysis
Si = sobol.analyze(problem, Y, print_to_console=True)

La sortie ressemble à ceci:

Parameter S1 S1_conf ST ST_conf
x1 0.307975 0.068782 0.560137 0.076688
x2 0.447767 0.057872 0.438722 0.042571
x3 -0.004255 0.052080 0.242845 0.025605

Parameter_1 Parameter_2 S2 S2_conf
x1 x2 0.012205 0.082296
x1 x3 0.251526 0.094662
x2 x3 -0.009954 0.070276

La sortie ressemble à ceci: L'indice de premier ordre, l'indice de second ordre et l'indice d'ordre total sont tracés comme suit pour une compréhension facile.

image.png image.png image.png

$ S_1 $ et $ S_2 $ ont de grandes valeurs, mais $ S_3 $ est presque nul. C'est parce que la fonction ʻIshigami` utilisée ici est définie comme suit, et $ S_3 $ de $ x_3 $ devient 0, mais il y a une interaction avec $ x_1 $, donc $ S_ {13} $ Sera une valeur positive.

f(x) = \sin(x_1) + a \sin^2(x_2) + b x_3^4 \sin(x_1)

Lorsque vous exécutez la méthode saltelli.sample, $ N (2D + 2) $ points d'échantillonnage sont générés, et vous devez trouver $ Y $ pour chacun. Si $ N $ est trop petit, la barre d'erreur sera grande, vous devez donc augmenter $ N $ jusqu'à ce que la barre d'erreur soit suffisamment petite. Généralement autour de N $ = 1000 $.

Cliquez ici pour le cahier Jupyter avec le code complet

Exécution par lots à l'aide d'OACIS

J'ai exécuté une fonction simple en python ici, mais en général, vous voudrez probablement utiliser un modèle mathématique plus lourd. Par exemple, si un calcul prend plusieurs minutes, il devient ingérable sur un PC et vous souhaitez l'exécuter sur une machine de cluster. Ici, nous allons introduire la procédure d'exécution par lots en utilisant oacis. Puisque oacis a une [API] Python (http://crest-cassia.github.io/oacis/en/api_python.html), vous pouvez l'utiliser pour soumettre des travaux et vous référer aux résultats.

Ici, à titre d'exemple, le script qui calcule la fonction Ishigami utilisée précédemment est enregistré en tant que simulateur. Après cela, si vous utilisez votre propre modèle mathématique, veuillez le lire comme il convient.

Tout d'abord, j'ai préparé un script pour calculer la fonction Ishigami comme suit. Il prend des paramètres comme arguments pour pouvoir être exécuté à partir d'oacis et écrit la sortie dans un fichier JSON appelé _output.json.

ishigami.py


import math,json

with open('_input.json') as f:
    x = json.load(f)
x1,x2,x3 = x['x1'],x['x2'],x['x3']
y = math.sin(x1) + 7.0*(math.sin(x2))**2 + 0.1*(x2**4)*math.sin(x1)
with open('_output.json', 'w') as f:
    json.dump({"y": y}, f)

Inscrivez-vous sur oacis avec les paramètres suivants.

Ensuite, préparez le script suivant.

from SALib.sample import saltelli
from SALib.analyze import sobol
from SALib.test_functions import Ishigami
import numpy as np

# import oacis module
import os,sys
sys.path.append( os.environ['OACIS_ROOT'] )
import oacis

# Define the model inputs
problem = {
    'num_vars': 3,
    'names': ['x1', 'x2', 'x3'],
    'bounds': [[-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359]]
}
# define problems
param_values = saltelli.sample(problem, 1000)
param_values

# find ishigami simulator and create ParameterSets and Runs
sim = oacis.Simulator.find_by_name('ishigami')
host = oacis.Host.find_by_name('localhost')
ps_array = []
for param in param_values:
    ps = sim.find_or_create_parameter_set( {"x1":param[0],"x2":param[1],"x3":param[2]} )
    ps.find_or_create_runs_upto( 1, submitted_to=host )
    ps_array.append(ps)

PS est fait avec ça. Après cela, oacis soumettra le travail sans autorisation, alors attendez qu'il soit terminé. (Ici, nous utilisons un modèle mathématique qui complète le calcul en un instant, mais comme oasis soumet chacun d'eux comme un travail, ce processus prend également environ une demi-journée. Cependant, un modèle mathématique qui prend plusieurs minutes ou plus. C'est très utile si vous souhaitez exécuter en parallèle sur un cluster etc., car une telle surcharge peut être ignorée)

image.png

Une fois terminé, si vous obtenez le résultat comme suit, le processus sera exactement le même.

# output values
Y = [ ps.average_result("y")[0] for ps in ps_array ]
Y = np.array(Y)
Y

Si = sobol.analyze(problem, Y)
Si

Cliquez ici pour le cahier jupyter de tous les codes utilisés cette fois

Recommended Posts

Analyse de sensibilité globale à l'aide des bibliothèques d'analyse de sensibilité Salib et Oacis
Analyse du réseau des acteurs vocaux (à l'aide de word2vec et networkx) (1/2)
Analyse du réseau des acteurs vocaux (utilisant word2vec et networkx) (2/2)