sobol全局灵敏性分析

1、Sobol全局灵敏性分析最近在研究全局敏感性分析方法中的Sobol方法，看了一些国内的论文，发现一个通病，就是公式 一挂就可以得出结果了，真心觉得这种论文很“恶心”，主要原因是自己看不太懂。直到在维基百科上面 找到了这种方法的详细解释，今天我们就根据网页上的步骤用一彳、例子来走-遍。1. 假设现在有一个函数：丫二sin(xj + 7 (sin(x2)2 + 0.1 x34 sin(xj，从公式中可以看到有x1、x2、x3三个自变量对应变量丫有影响。2. 然后一般会给这三个参数一个取值范围，这里假设三个自变量的取值范围都设为0,1。敏感性分析的目的就是求取这三个参数对于 丫值得贡献。当然我们这边

2、可能有人一下子就可以 分析出那个参数对于丫值影响最大，但是在解决实际问题时，这个函数一般都是未知，我们只能将其视 作一个黑盒子，只有输入和输出，这时我们对其进行敏感性分析就很有必要了。经过敏感性分析我们就能 找出对结果影响较大的参数。这样对于调整结果是很有帮助的。3. 接着上面的例子，首先我们得根据三个自变量的范围进行采样，这边采样的方法一般都是蒙特卡洛采样以及一系列基于蒙特卡洛采样的变种，这个例子中我们采用了 Sobol sequenee 具体的采样 原理在这就不说了，大家可以自行谷歌。为了方便讲解例子我们设置采样的样本数为4(N = 4)，自变量数目为3(D = 3)。我们按照上述网页的步

3、骤。4. 生成N*2D(即4行6列)的样本矩阵。这个就是我们Sobol sequenee做的事情。这边我们生成的矩 阵为：0.750.250.250.250.750.750.250.750.750.750.250.250.3750.3750.6250.8750.3750.1250.50.50.50.50.50.55、将矩阵的前D列设置为矩阵A、后D列设置为B列，在我们的例子中就是矩阵m的前3列设置为矩阵A,后3列设置为矩阵B o0.50.50.50.50.50.50.750.250.250.25c0.750.750.250.750.75B0.750.250.250.3750.3750.6250

4、.8750.3750.125构造NP的矩阵AB(i= 12Q厂即用矩阵B中的第i列替换矩阵A的第i列，以本体为例:0.50.50.50.50.50.50.50.50.50.250.250.250.750.750.250.750.250.75ABi0 7S0.750.75 AB20.250.25ABan 7S0.250.750.250.8750.3750.6250.3750.3750.6250.3750.3750.1256经过这三步我们构造了 A B、AB、ABAB这五个矩阵、这样我们就有(D + 2) * N (即20)组X儿、X3输入数据、因此我们将有20组丫值。将上述的数据带入函丫 = s

5、in(xj + 7(sin(X2)+ 0.1 X子sin(xj /这里详细的计算过程就不描述了。根据输入我们得出对应的丫值矩阵。2.091363877682.091363877682.091363877681.11036605855YaYb3.507651769323.507651769320.67596163481Yab11.11036605855a b 3.955626030581.310950363291.706651200031.718344245572.091363877682.091363877683.934324819331.1316672698Yaba0.68369300853

6、7 骨3.499920395591.310950363291.305370430237.根据一介影响指数公式：S Varxi Ex： (Y XJA/ar(Y)总效应指数曲皿昭叫師(丫)其中：t1 NiVarx.Ex： XJJ /侣办(f(AB)jf(A)Ex： .(Varxi(Y X： i)Var(Y) Var(YAYB)这里的不是普通意义上的相加而是构造了一个新的矩阵包括Ya和丫8即2.091363877681.110366058553.507651769321.31095036329Y2.091363877683.507651769321.110366058551.70665120003这

7、样我们就可以通过上述的公式求解x1的一介影响指数。下面是求解过程。8.Y的均值等于mea n(Y )=(2.+1.11036605855+3.50765176932+ 1.31095036329+2.+3.50765176932+ 1.11036605855+1.70665120003)/8 = 2.0545456218Var(Y) = 0.835332581542Varxi NjH(YABjYj)本例中 Varxi就等于(2(2. - 2.) + 3.50765176932* (0.67596163481 - 1.11036605855) +1.11036605855* (3.95562603058 3.50765176932) + 1.70665120003 * (1.71834424557 1.31095036329) )/4 = -0.62那么0.099075730.08276119315620.83533258154212N本例中 Ex“ 就等于(2. 2.)*2 + (1.11036605855 0.67596163481)9 + (3.50765176932 -3.95562603058卜2 + (1.31095036329 1.71834424557)八 2 )/(2 * 4)= 0.，那么St10.08310431220.069419739680.835