概要
Sobol
感度指数は、目的関数の分散分解に基づいてパラメータの重要度を定量化するグローバル感度分析手法です。
線形か非線形かを問わず、パラメータ間の交互作用も含めた影響度を [ 0 , 1 ] [0,1] [ 0 , 1 ]
の値として返します。
Tunny Dashboard では 2 種類の指数を提供します:
指数
記号
意味
一次指数 (First-Order)
S i S_i S i
パラメータ x i x_i x i の単独寄与率
全効果指数 (Total-Effect)
S T i ST_i S T i
x i x_i x i の単独 + 他パラメータとの交互作用を含む総寄与率
理論背景
ANOVA 分散分解
Sobol 分解は、モデル出力 Y = f ( X 1 , … , X p ) Y = f(X_1, \ldots, X_p) Y = f ( X 1 , … , X p )
の分散を各パラメータの寄与に分解します:
V a r ( Y ) = ∑ i V i + ∑ i < j V i j + ⋯ + V 1 … p \mathrm{Var}(Y) = \sum_i V_i + \sum_{i<j} V_{ij} + \cdots + V_{1\ldots p} Var ( Y ) = i ∑ V i + i < j ∑ V ij + ⋯ + V 1 … p
V i = V a r ( E [ Y ∣ X i ] ) V_i = \mathrm{Var}(\mathbb{E}[Y \mid X_i]) V i = Var ( E [ Y ∣ X i ]) : x i x_i x i のみが変化したときの Y Y Y
の期待値の分散
V i j = V a r ( E [ Y ∣ X i , X j ] ) − V i − V j V_{ij} = \mathrm{Var}(\mathbb{E}[Y \mid X_i, X_j]) - V_i - V_j V ij = Var ( E [ Y ∣ X i , X j ]) − V i − V j : x i x_i x i と
x j x_j x j の交互作用分散
一次指数
S i = V i V a r ( Y ) S_i = \frac{V_i}{\mathrm{Var}(Y)} S i = Var ( Y ) V i
x_i の変化だけで説明できる Y の分散の割合です。 S_i の合計は ≤
1(交互作用がなければ = 1)です。
全効果指数
S T i = E [ V a r ( Y ∣ X ∼ i ) ] V a r ( Y ) = 1 − V a r ( E [ Y ∣ X ∼ i ] ) V a r ( Y ) ST_i = \frac{\mathbb{E}[\mathrm{Var}(Y \mid X_{\sim i})]}{\mathrm{Var}(Y)} = 1 - \frac{\mathrm{Var}(\mathbb{E}[Y \mid X_{\sim i}])}{\mathrm{Var}(Y)} S T i = Var ( Y ) E [ Var ( Y ∣ X ∼ i )] = 1 − Var ( Y ) Var ( E [ Y ∣ X ∼ i ])
X_{~i} は「x_i 以外のすべてのパラメータ」です。 ST_i は、x_i の単独効果と x_i
が絡むすべての交互作用効果の合計です。
S T i ≥ S i ST_i \ge S_i S T i ≥ S i は常に成立します。 差 S T i − S i ST_i - S_i S T i − S i が大きいほど、 x i x_i x i
が他パラメータとの交互作用に強く関与していることを意味します。
有限サンプル推定の注意: Saltelli 一次指数推定量と Jansen
全効果推定量は独立に計算されるため、有限サンプルでは S ^ i > S T ^ i \hat{S}_i > \hat{ST}_i S ^ i > S T ^ i
が生じることがあります。 Tunny Dashboard では推定後に
S T ^ i ← max ( S T ^ i , S ^ i ) \hat{ST}_i \leftarrow \max(\hat{ST}_i, \hat{S}_i) S T ^ i ← max ( S T ^ i , S ^ i ) を適用してから [ 0 , 1 ] [0,1] [ 0 , 1 ]
にクランプし、理論的性質を出力値で保証しています。
Tunny Dashboard での計算方法
全体フロー
直接 Monte Carlo 積分を行うと n × p
倍のシミュレーション実行が必要になり、最適化問題では現実的ではありません。
代わりに二次 Ridge サロゲートモデルを構築し、そのモデル上で Saltelli
サンプリングを行います。
二次 Ridge サロゲートを構築します。
Saltelli 行列 A, B を ChaCha8 で生成します。
A, B, AB_i でサロゲートを評価します。
Jansen 推定量で S_i, ST_i を計算します。
Step 1: 二次 Ridge サロゲートモデル
特徴量構成
p 個の標準化済み線形入力から p(p+3)/2 次元の二次特徴量ベクトルを構築します:
ϕ ( x ) = [ x 1 , … , x p , x 1 2 , … , x p 2 , x 1 x 2 , … , x p − 1 x p ] \phi(x) = [x_1, \ldots, x_p,; x_1^2, \ldots, x_p^2,; x_1x_2, \ldots, x_{p-1}x_p] ϕ ( x ) = [ x 1 , … , x p , x 1 2 , … , x p 2 , x 1 x 2 , … , x p − 1 x p ]
例: p = 3 → 3 + 3 + 3 = 9 次元 (= 3×6/2)
したがって、サロゲートの特徴量次元は
d ϕ = p + p + p ( p − 1 ) 2 = p ( p + 3 ) 2 d_\phi = p + p + \frac{p(p-1)}{2} = \frac{p(p+3)}{2} d ϕ = p + p + 2 p ( p − 1 ) = 2 p ( p + 3 )
であり、学習サンプル数 n n n は少なくとも d ϕ d_\phi d ϕ と同程度、実務上は n ≫ d ϕ n \gg d_\phi n ≫ d ϕ が望ましいです。
サロゲートの学習
X X X の各列を Z スコア標準化
x j = x j − μ j σ j \tilde x_j = \frac{x_j-\mu_j}{\sigma_j} x j = σ j x j − μ j
二次特徴量 Φ \Phi Φ を構築
Φ \Phi Φ の各列を再度 Z スコア標準化
各目的関数に対して Ridge( α = 1.0 \alpha=1.0 α = 1.0 ) を適合
β k = ( Φ T Φ + I ) − 1 Φ T ( y k − y ˉ k ) \beta_k = (\Phi^T\Phi + I)^{-1}\Phi^T(y_k-\bar y_k) β k = ( Φ T Φ + I ) − 1 Φ T ( y k − y ˉ k )
サロゲートの評価
f ^ ( x ) = β T ϕ ( x ) + i n t e r c e p t \hat f(x) = \beta^T \tilde\phi(x) + \mathrm{intercept} f ^ ( x ) = β T ϕ ( x ) + intercept
ここで φ̃(x) は、二次特徴量を訓練データの統計量で標準化したものです。
Step 2: Saltelli サンプリング
疑似乱数生成: ChaCha8
ChaCha8 アルゴリズムに基づき、シードベースの決定論的な一様乱数([0, 1)
の範囲)を生成します。
初期シード: 0xDEAD_BEEF_1234_5678(固定・再現性あり)
行列 A, B の生成
各パラメータの値域 [lo_j, hi_j] は、実際のトライアルの min/max から取得します。
A [ s , j ] = l o j + u ( h i j − l o j ) , B [ s , j ] = l o j + u ( h i j − l o j ) A[s,j] = lo_j + u,(hi_j - lo_j),\qquad B[s,j] = lo_j + u,(hi_j - lo_j) A [ s , j ] = l o j + u ( h i j − l o j ) , B [ s , j ] = l o j + u ( h i j − l o j )
AB_i 行列
A の第 i 列のみを B の第 i 列で置換した行列です:
A B i [ s ] = [ A [ s , 0 ] , … , A [ s , i − 1 ] , B [ s , i ] , A [ s , i + 1 ] , … , A [ s , p − 1 ] ] AB_i[s] = [A[s,0], \ldots, A[s,i-1], B[s,i], A[s,i+1], \ldots, A[s,p-1]] A B i [ s ] = [ A [ s , 0 ] , … , A [ s , i − 1 ] , B [ s , i ] , A [ s , i + 1 ] , … , A [ s , p − 1 ]]
Step 3: サロゲート評価
f A [ k , s ] = s u r r o g a t e _ e v a l ( A [ s ] , o b j e c t i v e = k ) f_A[k,s] = \mathrm{surrogate_eval}(A[s],,\mathrm{objective}=k) f A [ k , s ] = surrogate_eval ( A [ s ] , objective = k )
f B [ k , s ] = s u r r o g a t e _ e v a l ( B [ s ] , o b j e c t i v e = k ) f_B[k,s] = \mathrm{surrogate_eval}(B[s],,\mathrm{objective}=k) f B [ k , s ] = surrogate_eval ( B [ s ] , objective = k )
f A B i [ k , s ] = s u r r o g a t e e v a l ( A B i [ s ] , o b j e c t i v e = k ) f{AB_i}[k,s] = \mathrm{surrogate_eval}(AB_i[s],,\mathrm{objective}=k) f A B i [ k , s ] = surrogate_eval ( A B i [ s ] , objective = k )
Step 4: Jansen 推定量
一次指数(Saltelli 2010)
S i ≈ ∑ s f B [ s ] ( f A B i [ s ] − f A [ s ] ) N V a r Y S_i \approx \frac{\sum_s f_B[s],(f_{AB_i}[s]-f_A[s])}{N,\mathrm{Var}_Y} S i ≈ N Var Y ∑ s f B [ s ] ( f A B i [ s ] − f A [ s ])
全効果指数(Jansen 1999)
S T i ≈ ∑ s ( f A [ s ] − f A B i [ s ] ) 2 2 N V a r Y ST_i \approx \frac{\sum_s (f_A[s]-f_{AB_i}[s])^2}{2N,\mathrm{Var}_Y} S T i ≈ 2 N Var Y ∑ s ( f A [ s ] − f A B i [ s ] ) 2
分散の推定
V a r Y = m e a n ( f A 2 ) − m e a n ( f A ) 2 \mathrm{Var}_Y = \mathrm{mean}(f_A^2) - \mathrm{mean}(f_A)^2 Var Y = mean ( f A 2 ) − mean ( f A ) 2
V a r Y ≈ 0 \mathrm{Var}_Y \approx 0 Var Y ≈ 0 の場合(サロゲートが定数に近い)は S i = S T i = 0 S_i = ST_i = 0 S i = S T i = 0
を返します。
クリッピング
推定誤差による範囲外値を防ぐため、最終値を [ 0 , 1 ] [0,1] [ 0 , 1 ] にクリップします:
S i = c l a m p ( S i r a w , 0 , 1 ) , S T i = c l a m p ( S T i r a w , 0 , 1 ) S_i = \mathrm{clamp}(S_i^{\mathrm{raw}}, 0, 1),\qquad ST_i = \mathrm{clamp}(ST_i^{\mathrm{raw}}, 0, 1) S i = clamp ( S i raw , 0 , 1 ) , S T i = clamp ( S T i raw , 0 , 1 )
複数目的関数への対応
Importance
チャートでは、各パラメータのスコアとして全目的関数に対する指数の平均を表示します:
d i s p l a y _ s c o r e ( p j ) = 1 m ∑ k S i [ j ] [ k ] \mathrm{display_score}(p_j)=\frac{1}{m}\sum_k S_i[j][k] display_score ( p j ) = m 1 k ∑ S i [ j ] [ k ]
d i s p l a y _ s c o r e ( p j ) = 1 m ∑ k S T i [ j ] [ k ] \mathrm{display_score}(p_j)=\frac{1}{m}\sum_k ST_i[j][k] display_score ( p j ) = m 1 k ∑ S T i [ j ] [ k ]
パラメータ設定
設定項目
値
Saltelli サンプル数
1024
Ridge 正則化強度 α
1.0
乱数シード
0xDEAD_BEEF_1234_5678
必要データ量
最低 2 トライアル(n ≥ 2)、p ≥ 1、目的関数 ≥ 1 が必要です
サロゲートの精度は n とともに向上します。目安として n ≥ 10 × p(p+3)/2
を推奨します
パラメータの扱いに関する前提
数値列はそのまま使用し、カテゴリ列は文字列ラベルを出現順の整数 ID(0.0, 1.0,
…)へラベル符号化して使用します。
カテゴリパラメータはラベル符号化済みの数値として扱われるため Sobol
計算の対象となりますが、整数 ID は順序や距離の情報を持たないため、Sobol
指数の解釈はカテゴリ変数に対して近似的なものとなります。
特性・限界
強み:
線形、非線形、交互作用をすべて扱えます
値域が [ 0 , 1 ] [0,1] [ 0 , 1 ] で複数パラメータ間の比較が容易です
ST_i - S_i で交互作用の強度を把握できます
サロゲートを使うため、トライアル数が限られていても計算可能です
弱み:
サロゲートモデル(二次 Ridge)の精度に結果が依存します
真の関数が高次非線形または強い不連続性を持つ場合は精度が落ちます
p が大きくなると特徴量数 p(p+3)/2
が増加し、サロゲートの学習精度が低下する可能性があります
乱数シードが固定のため再現性はありますが、N の選択で結果が変わります
参考文献
Saltelli, A. et al. (2010). Variance based sensitivity analysis of model
output. Design and estimator for the total sensitivity index. Computer
Physics Communications , 181(2), 259–270.
https://doi.org/10.1016/j.cpc.2009.09.018
Jansen, M. J. W. (1999). Analysis of variance designs for model output.
Computer Physics Communications , 117(1), 35–43.
https://doi.org/10.1016/S0010-4655(98)00154-4
Sobol, I. M. (1993). Sensitivity estimates for nonlinear mathematical models.
Mathematical Modelling and Computational Experiments , 1(4), 407–414.