【SALib】RBD-FASTを用いたグローバル感度解析

はじめに

修士論文執筆にあたりグローバル感度解析を扱う機会があり、それなりに勉強しましたので一旦記事にまとめます。
ご指摘、コメント等あればよろしくお願い致します。

本記事で紹介する内容

• ローカル感度解析・グローバル感度解析とは
• FAST（Forier Amlitude Sensitivity Test)について
• RBD-FASTについて
• SALibを用いた実装

感度解析

複数のパラメータを所有するモデルについて、そのパラメータが出力結果に対して与える影響を調べる手法が感度解析である。感度解析には多くの手法があり、大きく分けて「ローカル感度解析（Local Sensitivity Analysis）」、「グローバル感度解析（Global Sensitivity Analysis）」が挙げられる。

ローカル感度解析

ローカル感度解析は「one at a time（OAT）」と記載されることもあり、複数パラメータのうち、1つのパラメータを変動させ、その他のパラメータを固定した際に、変動させたパラメータが出力結果にどのような影響を与えるかを調べる手法である（pythonにて実装した例）。
この手法は計算コストが低い点、出力結果の変動の全てが変動させているパラメータに帰結するという点が利点として挙げられる一方で、パラメータの感度はあらかじめ設定された初期値に大きく依存する欠点も存在する。
Saltelli(2010)らによると、パラメータが10個であれば、ローカル感度解析にて探索できる範囲は全範囲の0.25%であり、特に非線形モデルである場合、ローカル感度解析は感度の結果に大きな不確実性を伴う。
また、ローカル感度解析の拡張であるMorris one-at-a-time法（グローバル感度解析の1つ）ではベースとなる初期値を複数選択し、アンサンブルな値を感度解析の結果とする（Morris, 1991;Wainwright, et.al.,2014）。

グローバル感度解析

グローバル感度解析はローカル感度解析と異なり、パラメータを固定することなく、変動する全てのパラメータの中で、それぞれのパラメータが出力結果に及ぼす全体に対する影響度を表す。グローバル感度解析は表に示す3種類に分類され(Tian, 2013)、本記事では特に分散法の中のFAST及び、FASTの派生形であるRBD-FASTについて取り扱う。

感度解析手法	代表手法	説明
回帰法	標準化回帰係数(SRC)、偏相関係数(PCC)	出力結果とパラメータについて回帰計算。 計算コストが低い。 非線形モデルの場合はランク変換して使用(SRCC, PCC)
スクリーニング法	Morris	出力結果の分散を減少させずに、入力パラメータの値を固定しながら感度解析。 最終的な感度はそれぞれのステップの平均値。 モデルに影響を及ぼす要因が少なく、影響の少ない因子が少ない場合に適用。 計算コストは低め。
分散法	FAST, Sobol, MARS, ACOSSO, SVM	各入力パラメータに対するモデル出力の分散を分解。主効果、相互作用効果の両方を考慮。 計算コストが高い。

分散法（ANOVA-based method）

グローバル感度解析における分散を基にした感度分析（ANOVA-based method）は、モデルの出力結果の分散に入力パラメータがどの程度影響しているかについて考える（例：あるモデルの出力結果の分散はその50%が降水量に、40%が日射量に、10%が降水量と日射量の相互作用によって構成されてる）。出力結果の分散への寄与度が入力パラメータの感度であり、本手法は非線形モデルにも適用可能である（Sobol, 2001）。
以降具体的な式を用い説明する（Goffart, et.al., 2015）。

あるモデル$y=f(x)$についてSobol分解を行う

y=f_0 + \sum_{i=1}^{m}f_i(x_i) + \sum_{i=1}^{m}\sum_{j=1}^{i-1}f_{ij}(x_i, x_j) + \dots + f_{1,\ldots,m}(x_1, x_2, \ldots, x_m)

このとき
$f_0$:モデル出力結果の平均値
$f_i(x_i)$:モデル出力の$x_i$にのみ起因する部分
$f_{ij}(x_i, x_j)$: モデル出力の$x_i$, $x_j$の両方に影響される部分
$f_{1,\ldots,m}(x_1, x_2, \ldots, x_m)$: 全パラメータによって影響を受ける部分

この時、分散$V$は以下のようにあらわすことが出来る。

V=\sum_{i=1}^{m}V_i + \sum_{i=1}^{m}\sum_{j=1}^{i-1}V_{ij} + \dots + V_{1,\ldots,m}

さらに両辺を$V$で除し、相対感度$S$として表す。

1=\sum_{i=1}^{m}S_i + \sum_{i=1}^{m}\sum_{j=1}^{i-1}S_{ij} + \dots + S_{1,\ldots,m}

このとき、
$V$: 出力結果の分散
$V_i$: $x_i$の変動のみが出力結果に及ぼす分散
$V_{ij}$: $x_i$, $x_j$の両方の変動によってのみ出力結果に及ぼす分散
$V_{1,\ldots,m}$: 全パラメータによってのみ出力結果に及ぼす分散

$S$についても同様であり、
特に$S_i$は$x_i$の主感度であり$x_i$のみが及ぼす1次の感度指数であり、2次からm次までの感度指数がモデル中の入力の相互作用によって説明される部分である。したがってパラメータ間に相互作用が見られない場合以下の式を満たす。

\sum_{i=1}^{m}S_i=1

分散に基づくグローバル感度解析は、上記の考え方に基づき感度を求める。そして、分散を求める手法としてフーリエ変換を用いた手法がFASTである。

Forier Amlitude Sensitivity Test（FAST）

フーリエ振動感度テストは、感度の高いパラメータの周波数が出力結果の周波数にも強く表現されるという前提に基づく感度解析手法である。
図1（Kraft and Rinderknecht(2018)より引用）に示すように、感度解析対象パラメータのサンプリングはランダムではなく、それぞれのパラメータに固有の周波数に従って行い、出力結果の周波数を解析し、その周波数及び整数倍周波数に対応した高次ハーモニクスのパワースペクトルをそのパラメータの感度とする手法である。

以降具体的な式変形を用いて説明していく（劉, 2010; Saltelli,et.al., 2012）。

原点回りの$r$次モーメント$\langle y^{(r)}\rangle$について以下の式にて定義できる（$x$は単位超立方体に一様に分布していると仮定）。

\langle y^{(r)}\rangle=\idotsint_{\Omega_n}f^r(x_1, x_2,\cdots,x_n) dx_1,dx_2,\cdots,dx_n

$\Omega_n$は、すべての入力変数が構成する$n$次元空間であり、$y$の分散$V_y$について以下の式が成り立つ。

V_y=\langle y^{(2)}\rangle-(\langle y^{(1)}\rangle)^2

次に、以下のように$x_i$を変換する。

x_i(s)=G_i(\sin(\omega_is))

ここで$\omega_i$はそれぞれの$x_i$に割り当てられた固有周波数、$s$は実数であり$s$を変動させることで、入力変数$x_1, x_2,\cdots,x_n$が$\Omega_n$空間で同時に変動し、この変換によって、$n$次元の入力空間を1次元の空間として扱うことができ、計算が簡便となる。
$x_i$は$\omega_i$に従って振動し、どのようなモデル$f$であってもその結果はそれぞれ異なる周波数の組み合わせとして現れる。もし$i$番目の入力が強い影響を結果に及ぼしていた場合、$\omega_i$における$y$の振動も強く表れる。また、それぞれのパラメータは高次ハーモニクスにて干渉する。
$G_i$はエルゴード性を満たす探索関数として設定され、例として$G_i=\frac{1}{\pi}\arcsin[\sin(\omega_is)]+0.5$が挙げられる。

変換後の式についてエルゴード定理より$s$が無限区間で変動するとき

\begin{align}
\bar{y^{(r)}} &= \lim_{x \to \infty} \frac{1}{2T} \int^{T}_{-T}f^r(x_1(s), x_2(s),\ldots,x_n(s))ds \\
&= \langle y^{(r)} \rangle
\end{align}

さらに$\omega_i\in \mathbb{Z}$を満たすとすると$\sin(\omega_is)=\sin(\omega_is+2\omega_i\pi)$を満たすため、以下のように式変形できる。

\begin{align}
\bar{y^{(r)}} &= \lim_{x \to \infty} \frac{1}{2T} \int^{T}_{-T}f^r(x_1(s), x_2(s),\ldots,x_n(s))ds \\
&= \frac{1}{2\pi} \int^{\pi}_{-\pi}f^r(x_1(s), x_2(s),\ldots,x_n(s))ds \\
\end{align}

したがって

\begin{align}
V_y &= \langle y^{(2)}\rangle-(\langle y^{(1)}\rangle)^2 \\
&=\bar{y^{(2)}}-(\bar{y^{(1)}})^2 \\
&=\frac{1}{2\pi} \int^{\pi}_{-\pi}f^2(x_1(s), x_2(s),\ldots,x_n(s))ds-[\frac{1}{2\pi} \int^{\pi}_{-\pi}f(x_1(s), x_2(s),\ldots,x_n(s))ds]^2
\end{align}

ここで$f(s)$についてフーリエ変換を行うと

\begin{align}
&y=f(x_1(s), x_2(s),\ldots,x_n(s))ds=\sum_{j\in \mathbb{Z}}[A_j\cos(js)+B_j\sin(js)] \\
&A_j=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x_1(s), x_2(s),\ldots,x_n(s))ds\cos(js)ds \\
&B_j=\frac{1}{2\pi}\int_{-\pi}^{\pi}f(x_1(s), x_2(s),\ldots,x_n(s))ds\sin(js)ds
\end{align}

したがって$\bar{y^{(2)}}$は

\begin{align}
\bar{y^{(2)}}&=\frac{1}{2\pi}\int_{-\pi}^{\pi}f^2(x_1(s), x_2(s),\ldots,x_n(s))ds \\
&=\frac{1}{2\pi}\int_{-\pi}^{\pi}[\sum_{j\in \mathbb{Z}}[A_j\cos(js)+B_j\sin(js)] ]^2ds \\
&=\sum_{j\in \mathbb{Z}}(A_j^2+B_j^2)
\end{align}

よって$A_{-j}=A_j$,$B_{-j}=B_j$および以下の式が成り立つ。

\begin{align}
A_0&=\frac{1}{2\pi}\int^{\pi}_{-\pi}f(x_1(s), x_2(s),\ldots,x_n(s))\cos0ds \\
&=\frac{1}{2\pi}\int^{\pi}_{-\pi}f(x_1(s), x_2(s),\ldots,x_n(s))ds \\
&=\langle y^{(1)}\rangle
\end{align}

そのため

\begin{align}
V_y &= \langle y^{(2)}\rangle-(\langle y^{(1)}\rangle)^2 \\
&=\sum_{j\in \mathbb{Z}}(A_j^2+B_j^2)-A_0^2 \\
&=2\sum_{j\in \mathbb{Z}^+}(A_j^2+B_j^2)+A_0^2-A_0^2 \\
&=2\sum_{j\in \mathbb{Z}^+}(A_j^2+B_j^2)
\end{align}

ゆえに$x_i$の$y$への主効果$S_i$は以下の式にて求められる。

\begin{align}
S_i &= \frac{V_i}{V} \\
V_i &= 2\sum_{p=1}^{+\infty}(A_{p\omega_i}^2 + B_{p\omega_i}^2) \\
V &= 2\sum_{j\in \mathbb{Z}^+}(A_j^2 + B_j^2)
\end{align}

したがって$x_i$の1次感度は割り当てた特性周波数($\omega_i$)及び高次ハーモニクス($2\omega_i, 3\omega_i \cdots$)のパワースペクトルが全体のパワースペクトルに占める割合として計算され、任意の相互作用はそれ以外の周波数におけるパワースペクトルに含まれる。

Random Balanced Design・RBD-FAST

先に説明したFASTは複数の周波数を用いることから、サンプル数$N$に制限があり、干渉を起こさない周波数についての最大次数$M$、固有周波数$\omega_i$から以下式のように定まる。

N\geq2M max(\omega_i)

この欠点を克服する手法として、Random Balanced Design(RBD)、およびRBD-FASTがあり(Tissot and Prieur, 2012)、これらの手法はパラメータ数の増加に伴うサンプル数の制限が無い（要はFASTに比べて少ないサンプル数で感度解析ができる）。RBD及びRBD-FASTでは、FASTと異なり、それぞれのサンプルに対してそれぞれの周波数を割り当てず、RBDでは全てのパラメータについて同様の周波数を割り当て、RBD-FASTでは、入力パラメータをいくつかのグループに分割し、それぞれのグループについて周波数を割り当てる。RBD-FASTとRBDの違いは、RBD-FASTでは任意の次数の感度解析が、RBDでは一次の感度解析のみが得られるという点である(Tissot and Prieur, 2012)。従って、RBDはRBD-FASTに包含されていると考えられる（1つの周波数のみを割り当てるように分割した場合）。
また、SALibにおけるRBD-FASTでは1次の感度解析のみを行うEASI RBD-FASTを使用しており、以降EASI RBD-FASTについて解説する（Goffart and Woloszyn. 2021）。

EASI RBD-FAST

EASI RBD-FASTはRBD-FASTについてサンプル整形にEASIを使用した感度解析であり、1次の感度しか解析できない一方で、既存の出力結果についても適用できる利点を持つ。図2(Goffart and Woloszyn(2021)より引用, 改変)に概要を示す。

パラメータサンプリング

それぞれのパラメータを特定の周波数に従ってサンプリングし、それをランダム順に並べ替え、それらによって作成したパラメータセットでモデルを実行し出力結果を得る。周波数は一般的に最も簡単な$\omega_i=1$を設定する。
従って、特定の周波数に並べ替えられるような、パラメータセットが作成できていればよいので、モンテカルロ法、準モンテカルロ法、ラテン超方格法等でランダムサンプリングを行なえばよい（FASTと違いパラメータサンプリングに強い制約はない）。従って、既に出力済みのシミュレーション結果を感度解析に使用することができる。

サンプル数

Goffart and Woloszyn(2021) はサンプル数は数百回程度で十分と述べており、適切なサンプル数を推定する手法として、モデル出力結果の平均値と分散についてシミュレーション回数毎にプロットを行い、平均値や分散の変動が安定したシミュレーション回数を適切なサンプル数とする手法を挙げている。

感度解析

パラメータそれぞれについてのランダムサンプルを特定の周波数の三角波形に変換し、それに伴って並べ替えられる出力結果についてフーリエ感度解析を行う（従ってフーリエ変換に伴う感度解析をパラメータの種類ぶん行う）。$\omega_i=1$ではランダムサンプリング後のパラメータを昇順に並べ、それらについて偶数番目を抜き出して降順にして、昇順の奇数番目の配列の後ろに結合することで三角波形を作成する。この整形を、Effective Algorithm for variance-based Sensitivity Indices (EASI) と呼ぶ。
それぞれの感度解析では得られたパワースペクトルがパラメータの1次感度であり、注目していないパラメータ（三角波形になっていない）はホワイトノイズとして扱われる。加えて感度推定精度向上のために、計算後感度にホワイトノイズの補正を考慮した式が以下である（Tissot and Prieur. 2012）。なお、バイアス自体はサンプルサイズが大きく感度が高いほど影響が小さくなる。

\begin{align}
S_{revised} &= S_1 - \frac{\lambda}{1-\lambda} \\
\lambda&=\frac{2M}{N} 
\end{align}

$M$は合計する高次ハーモニクスの次数、$N$はサンプルサイズである

SALibによるRBD-FAST

SALibに実装されたRBD-FASTを実際に使用してみる
このコードではラテン超方格にてサンプリング後、モデルを実行し出力値とともに感度解析を行っている。標準偏差(conf)はブーストトラップサンプリングによって計算される。サンプル数が足りていないと標準偏差が大きくなる。

from SALib.analyze import rbd_fast
from SALib.sample import latin
from SALib.test_functions import Ishigami

problem = {
    'num_vars': 3,
    'names': ['x1', 'x2', 'x3'],
    'bounds': [[-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359],
               [-3.14159265359, 3.14159265359]]
}


X = latin.sample(problem, 1000)
Y = Ishigami.evaluate(X)
Si = rbd_fast.analyze(problem, X, Y, print_to_console=True)

結果は以下。マイナスの値が存在するのはバイアス補正が原因?

          S1   S1_conf
x1  0.296304  0.058937
x2  0.448655  0.071417
x3 -0.003163  0.035991

---

勿論先に説明した通り、既に実行済みのモデル計算結果を利用し、直接rbd_fast.analyzeを利用してもよい。以下はcsvからデータを読み取り感度解析を行った例。

from SALib.analyze import rbd_fast
import pandas as pd
import numpy as np
problem = {
    'num_vars': 3,
    'names': ['x1', 'x2', 'x3'],
    'bounds': [[0,1],
               [0,1],
               [0,1]]
}
df = pd.read_csv('test_data.csv')
X = df.iloc[:,:-2].values
Y = df.iloc[:,-1].values
Si = rbd_fast.analyze(problem, X, Y, print_to_console=True)

読み込むcsvの形式

x1	x2	x3	y
0.54	0.71	0.60	2.55
0.42	0.64	0.43	2.01
0.96	0.38	0.79	3.17
0.56	0.92	0.07	2.48
...	...	...	...

おわりに

今回はデータが使いまわせるFASTという点に魅力を感じ、RBD-FASTについて掘り下げてみました。1次の感度しか計算できないので2次以上の感度となると別の感度解析を試す必要がありますが、SALibで簡単に使えますので、かなりとっつきやすいです。

参考文献