LoginSignup
23
13

More than 3 years have passed since last update.

@yohm(Yohsuke Murase)

感度解析ライブラリsalibとoacisを使った大域的感度解析

Posted at

感度解析(Sensitivity analysis)とは？

感度解析（Sensitivity analysis）とは、数理モデルやシステムの出力の不確実さ(uncertainty)が、どの入力の不確実さに由来するかを定量的に評価するための手法。
https://en.wikipedia.org/wiki/Sensitivity_analysis

出力$Y$、入力$\bf X$とした時に数理モデル以下のようにあらわせるとする。ある入力$X_i$に不確実性があると$Y$にも不確実性が現れるが、$Y$の不確実性がどの$X_i$の不確実性にどれくらい依存しているかを定量的に評価する。

Y = f({\bf X})

以下のような目的に使用される。

  • モデルに不確実性があるとき（たとえば直接観測できない不明確な量を仮定した場合など）に、得られた結果がどれくらい頑健かをテストする
  • モデルの入力と出力の間の関係性を理解する
  • 感度の大きい入力を特定し、その入力に注力して確実性を上げることで計算結果の確実性を向上させる
  • 感度の小さい入力を固定することによるモデルの単純化
  • 多くのパラメータを持つモデルのパラメータをキャリブレートするときに、どの入力に注力してキャリブレートするべきかを決めるのに使う

感度解析の手法

感度解析にはいくつかの手法がある。計算コストや目的、数理モデルの特徴によって使い分ける。

もっとも単純な手法としてOne-at-a-time (OAT)がある。名前の通り、一つのパラメータだけを動かしてみて結果がどれくらい変わるかを偏微分係数や線形回帰などの手法で調べる手法である。とても単純でわかりやすいし計算コストも小さいが、動かすパラメータ以外のパラメータをベースラインの値に固定してしまうので、入力空間をフルに探索できない。

また、偏微分係数で評価しようとすると、その特定のベースラインでの値を局所的にしか評価することができない。$y$が$x_i$に対して非単調に変化する場合や入力変数間の相互作用を評価できない

global sensitivity analysis (Sobol method)

これらの問題を解決する大域的な感度解析(global sensitivity analysis)を行う手法としてVariance-based sensitivity analysis (分散ベースの感度解析）がある。よくSobol methodと呼ばれたりもする。
https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis

この手法では入力、出力を確率変数として考える。出力の分散を、ある一つの入力（もしくは複数の入力）による寄与に分解する。例えば、出力の分散の70%が$X_1$のゆらぎ由来で、20%が$X_2$、残り10%が$X_1$と$X_2$の相互作用による寄与であるということが言える。

この手法は(i)大域的、(ii)非線形な反応も扱うことができ、(iii) 入力パラメータの相互作用（複数のパラメータを同時に変化させた時に現れる効果）を測定できるという手法になっている。
一方で統計的に処理するので、入力空間を十分にサンプリングする必要があり、一般に様々な入力での計算が多数必要になり計算コストが重い。

具体的にどのような量を計算するか。今モデルを$d$次元の変数$\bf X = {X_1,X_2,\dots,X_d}$とする。
また全ての変数は$X_i \in [0,1]$と規格化され均一に分布しているものとする。出力$Y$は一般的に以下のようにかくことができる。

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

ここで$f_i$は$X_i$の関数、$f_{ij}$は$X_i$,$X_j$の関数である。

$Y$の分散を計算すると

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

と書くことができ、分散が各$i$に対して依存する項の和として書くことができる。ここで$V_i$は以下の式で定義される。

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

この式は、まず$X_i$を特定の値に固定して$i$以外の変数に対して積分して特定の$X_i$に対する$Y$の期待値を計算し、その量の$X_i$を変動させたときの分散を計算することを意味している。

二次の項$V_{ij}$も同様に定義されて以下のようになる。

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}

さてパラメータ$i$の感度を評価する指標として主にfirst-order index $S_i$, total-effect index $S_{Ti}$の２種類がある。それぞれ以下のように定義される。

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$は$i$「のみ」を変化させたときの$Y$の分散を評価する手法である。$i$単独で変動させた時に平均して$S_i$ほど揺らぐ。
一方、Total-effect indexはすべてのパラメータを変動させたときの$i$の寄与分で、他のパラメータとの相互作用の寄与も含めたもの。$S_{Ti}$が小さければ、その値は固定しておいても結果に影響を与えないのでモデルの単純化などに利用できる。

これらの指標を計算するために、パラメータ空間をモンテカルロサンプリングする。
ランダムサンプリングをしてもよいが、より効率的に計算するための工夫として、Quasi-Monte Carlo methodを使う場合が多く、実際以下で説明するライブラリでもそうしている。通常、数千x(パラメータ数)程度のサンプリングが必要になる。

SALibの使い方

global sensitivity analysisを行うライブラリの一つにSALibがある。
ここでは使い方を説明する。
このライブラリにはSobol methodをはじめとしていくつかの手法が実装されているが、ここではSobol methodのみを解説する。

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

インストール方法はpip install SALib
前提条件としてnumpy,scipy,matplotlibが必要。

SALibは以下の４ステップで行う。

  1. モデルインプットの定義
  2. サンプルの生成
  3. 感度解析したい数理モデルの実行
  4. 感度解析の指標($S_i$や$S_{Ti}$)を計算

もっとも簡単なサンプルコードは以下のようになる。
ここではIshigamiというテスト用の関数を数理モデルとして利用している。

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)

出力は以下のようになる。

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

出力は以下のようになる。わかりやすくなるようにfirst-order index, second-order index, total-order indexをプロットしてみると以下のようになる。

image.png
image.png
image.png

$S_1$,$S_2$は大きな値が出ているが、$S_3$はほぼゼロである。
これは、ここで使っているIshigami関数は以下のように定義されており、$x_3$の$S_3$は0になるが、$x_1$との相互作用があるので$S_{13}$は正の値になる。

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

saltelli.sampleメソッドを実行すると$N(2D+2)$個のサンプル点が生成され、それぞれについて$Y$を求める必要がある。$N$が小さすぎるとエラーバーが大きくなるので、エラーバーが十分に小さくなるまで$N$を増やす必要がある。一般的には$N=1000$くらい。

全コードを含むjupyter notebookはこちら

OACISを使ったバッチ実行

ここでは単純な関数をpythonで実行したが、一般的にはより計算の重い数理モデルを使うことになるだろう。例えば一つの計算に数分かかるとすると一つのPCでは手に負えなくなり、クラスターマシンなどで実行したくなる。
ここではoacisを使ってバッチ実行する手順を紹介する。oacisにはPythonのAPIがあるので、それを使ってジョブ投入や結果の参照を行う。

ここでは例として先ほど使ったIshigami関数を計算するスクリプトをsimulatorとして登録する。以降、自分の数理モデルを使う場合は適宜読み替えて欲しい。

まずIshigami関数を計算するスクリプトを以下のように用意した。oacisから実行できるようにパラメータを引数で受け取り、出力を_output.jsonという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)

oacisには以下のような設定で登録する。

  • Name: Ishigami
  • Definition of parameters: [(x1,Float), (x2,Float), (x3,Float)]
  • command: "python ~/path/to/ishigami.py"

続いて、以下のようなスクリプトを用意する。

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が作られる。あとは勝手にoacisがジョブを投入してくれるので完了するまで待つ。
（ここでは一瞬で計算が終わる数理モデルを使っているが、oacisは一つ一つをジョブとして投入するのでこの処理にも半日ほどの時間がかかる。しかし、数分以上の時間のかかる数理モデルをクラスターなどで並列に実行したい場合はそのようなオーバーヘッドは無視できるので、非常に有用である）

image.png

完了後、以下のように結果を取得すれば、以後は全く同じ処理になる。

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

Si = sobol.analyze(problem, Y)
Si

今回利用した全コードのjupyter notebookはこちら

23
13
0

Register as a new user and use Qiita more conveniently

  1. You get articles that match your needs
  2. You can efficiently read back useful information
  3. You can use dark theme
What you can do with signing up
Sign upLogin
23
13