[PYTHON] Globale Sensitivitätsanalyse unter Verwendung der Sensitivitätsanalysebibliotheken Salib und Oacis

Was ist Sensitivitätsanalyse?

Die Sensitivitätsanalyse ist eine Methode zur quantitativen Bewertung, welche Eingangsunsicherheit aus der Unsicherheit der Ausgabe eines mathematischen Modells oder Systems abgeleitet wird. https://en.wikipedia.org/wiki/Sensitivity_analysis

Wenn die Ausgabe $ Y $ und die Eingabe $ \ bf X $ ist, wird das mathematische Modell wie folgt ausgedrückt. Wenn eine Eingabe $ X_i $ unsicher ist, erscheint die Unsicherheit auch in $ Y $, aber quantitativ hängt die Unsicherheit in $ Y $ davon ab, welche $ X_i $ unsicher ist. bewerten.

Y = f({\bf X})

Es wird für die folgenden Zwecke verwendet.

Sensitivitätsanalysemethode

Es gibt verschiedene Methoden zur Sensitivitätsanalyse. Verwenden Sie es ordnungsgemäß gemäß den Berechnungskosten, dem Zweck und den Merkmalen des mathematischen Modells.

One-at-a-Time (OAT) ist die einfachste Methode. Wie der Name schon sagt, können Sie überprüfen, um wie viel sich das Ergebnis ändert, indem Sie nur einen Parameter mit einer Methode wie dem partiellen Differentialkoeffizienten oder der linearen Regression verschieben. Es ist sehr einfach und leicht zu verstehen und die Berechnungskosten sind gering. Da jedoch andere Parameter als die Bewegungsparameter auf die Basislinienwerte festgelegt sind, kann der Eingaberaum nicht vollständig durchsucht werden.

Wenn versucht wird, anhand des partiellen Differentialkoeffizienten zu bewerten, kann der Wert an der spezifischen Basislinie nur lokal bewertet werden. Die Interaktion zwischen Eingabevariablen kann nicht ausgewertet werden oder wenn sich $ y $ in Bezug auf $ x_i $ nicht monoton ändert

global sensitivity analysis (Sobol method)

Die varianzbasierte Sensitivitätsanalyse ist eine Methode zur Durchführung einer globalen Sensitivitätsanalyse zur Lösung dieser Probleme. Es wird oft als Sobol-Methode bezeichnet. https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis

Bei dieser Methode werden Eingabe und Ausgabe als stochastische Variablen betrachtet. Teilen Sie die Verteilung der Ausgabe in Beiträge von einer Eingabe (oder mehreren Eingaben) auf. Zum Beispiel kann gesagt werden, dass 70% der Ausgabeverteilung aus der Fluktuation von $ X_1 $ stammen, 20% aus $ X_2 $ und die restlichen 10% aus der Interaktion zwischen $ X_1 $ und $ X_2 $.

Diese Methode kann (i) globale und (ii) nichtlineare Reaktionen verarbeiten und (iii) die Wechselwirkung von Eingabeparametern messen (Effekte, die auftreten, wenn mehrere Parameter gleichzeitig geändert werden). Es gibt. Andererseits ist es, da es statistisch verarbeitet wird, notwendig, den Eingaberaum ausreichend abzutasten, und im Allgemeinen ist eine große Anzahl von Berechnungen mit verschiedenen Eingaben erforderlich, und die Berechnungskosten sind hoch.

Welche Art von Betrag sollte speziell berechnet werden? Nun sei das Modell eine $ d $ -Dimensionsvariable $ \ bf X = \ {X_1, X_2, \ dots, X_d } $. Es wird angenommen, dass alle Variablen als $ X_i \ in [0,1] $ standardisiert und gleichmäßig verteilt sind. Die Ausgabe $ Y $ kann im Allgemeinen wie folgt geschrieben werden:

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})

Wobei $ f_i $ eine Funktion von $ X_i $ ist und $ f_ {ij} $ eine Funktion von $ X_i $ ist, $ X_j $.

Wenn Sie die Verteilung von $ Y $ berechnen

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

Kann als die Summe der Begriffe geschrieben werden, deren Varianz von jedem $ i $ abhängt. Hier wird $ V_i $ durch die folgende Formel definiert.

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

Diese Formel fixiert zuerst $ X_i $ auf einen bestimmten Wert und integriert ihn in andere Variablen als $ i $, um den erwarteten Wert von $ Y $ für einen bestimmten $ X_i $ und diesen Betrag von $ X_i $ zu berechnen Es bedeutet, die Varianz zu berechnen, wenn die geändert wird.

Der quadratische Term $ V_ {ij} $ ist ebenfalls wie folgt definiert.

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}

Es gibt hauptsächlich zwei Arten von Indizes zur Bewertung der Empfindlichkeit des Parameters $ i $: den Index erster Ordnung $ S_i $ und den Gesamteffektindex $ S_ {Ti} $. Jedes ist wie folgt definiert.

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 $ ist eine Methode zur Bewertung der Verteilung von $ Y $, wenn $ i $ "nur" geändert wird. Wenn $ i $ alleine geändert wird, schwankt es im Durchschnitt um ungefähr $ S_i $. Andererseits ist der Total-Effect-Index der Beitrag von $ i $, wenn alle Parameter geändert werden, einschließlich des Beitrags der Interaktion mit anderen Parametern. Wenn $ S_ {Ti} $ klein ist, wirkt sich der Wert nicht auf das Ergebnis aus, selbst wenn er fest ist, sodass er zur Vereinfachung des Modells verwendet werden kann.

Zur Berechnung dieser Indikatoren wird der Parameterraum mit Monte Carlo abgetastet. Zufallsstichproben können verwendet werden, aber die Quasi-Monte-Carlo-Methode wird häufig als Hilfsmittel für eine effizientere Berechnung verwendet. Tatsächlich tun dies die unten beschriebenen Bibliotheken. Normalerweise ist eine Abtastung von mehreren tausend x (Anzahl der Parameter) erforderlich.

Wie benutzt man SALib?

Eine der Bibliotheken, die eine globale Sensitivitätsanalyse durchführen, ist SALib. Hier wird erklärt, wie man es benutzt. In dieser Bibliothek sind mehrere Methoden implementiert, einschließlich der Sobol-Methode. Hier wird jedoch nur die Sobol-Methode erläutert.

Referenz: https://www.youtube.com/watch?v=gkR_lz5OptU&feature=youtu.be

Die Installationsmethode ist "pip install SALib" Benötigt numpy, scipy, matplotlib als Voraussetzungen.

SALib wird in den folgenden 4 Schritten ausgeführt.

  1. Definition der Modelleingabe
  2. Probengenerierung
  3. Ausführung des mathematischen Modells, für das eine Sensitivitätsanalyse gewünscht wird
  4. Berechnen Sie die Sensitivitätsanalyseindikatoren ($ S_i $ und $ S_ {Ti} $).

Der einfachste Beispielcode lautet wie folgt. Hier wird eine Testfunktion namens "Ishigami" als mathematisches Modell verwendet.

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)

Die Ausgabe sieht folgendermaßen aus:

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

Die Ausgabe sieht folgendermaßen aus: Der Index erster Ordnung, der Index zweiter Ordnung und der Index totaler Ordnung sind zum leichteren Verständnis wie folgt dargestellt.

image.png image.png image.png

$ S_1 $ und $ S_2 $ haben große Werte, aber $ S_3 $ ist fast Null. Dies liegt daran, dass die hier verwendete Ishigami-Funktion wie folgt definiert ist und $ S_3 $ von $ x_3 $ zu 0 wird, aber es gibt eine Interaktion mit $ x_1 $, also $ S_ {13} $ Wird ein positiver Wert sein.

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

Wenn Sie die Methode saltelli.sample ausführen, werden $ N (2D + 2) $ Beispielpunkte generiert, und Sie müssen für jeden $ Y $ finden. Wenn $ N $ zu klein ist, ist die Fehlerleiste groß. Sie müssen also $ N $ erhöhen, bis die Fehlerleiste klein genug ist. Im Allgemeinen um $ N = 1000 $.

Klicken Sie hier für ein Jupyter-Notizbuch mit vollständigem Code

Stapelausführung mit OACIS

Ich habe hier eine einfache Funktion in Python ausgeführt, aber im Allgemeinen möchten Sie wahrscheinlich ein schwereres mathematisches Modell verwenden. Wenn eine Berechnung beispielsweise mehrere Minuten dauert, kann sie auf einem PC nicht mehr verwaltet werden und Sie möchten sie auf einem Cluster-Computer ausführen. Hier werden wir das Verfahren für die Stapelausführung mit oacis vorstellen. Da oacis über eine Python [API] verfügt (http://crest-cassia.github.io/oacis/en/api_python.html), können Sie damit Jobs senden und auf die Ergebnisse verweisen.

Hier wird beispielsweise das Skript, das die zuvor verwendete Ishigami-Funktion berechnet, als Simulator registriert. Wenn Sie danach Ihr eigenes mathematisches Modell verwenden, lesen Sie es bitte entsprechend.

Zuerst habe ich ein Skript vorbereitet, um die Ishigami-Funktion wie folgt zu berechnen. Es verwendet Parameter als Argumente, damit es von oacis ausgeführt werden kann, und schreibt die Ausgabe in eine JSON-Datei mit dem Namen "_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)

Registrieren Sie sich bei oacis mit den folgenden Einstellungen.

Bereiten Sie als Nächstes das folgende Skript vor.

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 ist damit gemacht. Danach sendet oacis den Job ohne Erlaubnis. Warten Sie also, bis er abgeschlossen ist. (Hier verwenden wir ein mathematisches Modell, das die Berechnung sofort abschließt. Da Oasis jedoch jedes einzelne als Job einreicht, dauert dieser Vorgang auch etwa einen halben Tag. Bei einem mathematischen Modell dauert es jedoch mehrere Minuten oder länger. Dies ist sehr nützlich, wenn Sie parallel auf einem Cluster usw. ausgeführt werden möchten, da dieser Overhead ignoriert werden kann.

image.png

Wenn Sie nach Abschluss das folgende Ergebnis erhalten, ist der Vorgang genau der gleiche.

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

Si = sobol.analyze(problem, Y)
Si

Klicken Sie hier, um das Jupyter-Notizbuch mit allen diesmal verwendeten Codes anzuzeigen

Recommended Posts

Globale Sensitivitätsanalyse unter Verwendung der Sensitivitätsanalysebibliotheken Salib und Oacis
Netzwerkanalyse von Sprachschauspielern (mit word2vec und networkx) (1/2)
Netzwerkanalyse von Sprachschauspielern (mit word2vec und networkx) (2/2)