sobol全局灵敏性分析

Sobol全局灵敏性分析最近在研究全局敏感性分析方法中的Sobol方法，看了一些国内的论文，发现一个通病，就是公式一挂就可以得出结果了，真心觉得这种论文很“恶心”，主要原因是自己看不太懂。直到在维基百科上面找到了这种方法的详细解释，今天我们就根据网页上的步骤用一个例子来走一遍。1.假设现在有一个函数:，从公式中可以看到有x1、x2、x3三个自变量对应变量Y有影响。2.然后一般会给这三个参数一个取值范围，这里假设三个自变量的取值范围都设为[0,1]。敏感性分析的目的就是求取这三个参数对于Y值得贡献。当然我们这边可能有人一下子就可以分析出那个参数对于Y值影响最大，但是在解决实际问题时，这个函数一般都是未知，我们只能将其视作一个黑盒子，只有输入和输出，这时我们对其进行敏感性分析就很有必要了。经过敏感性分析我们就能找出对结果影响较大的参数。这样对于调整结果是很有帮助的。3.接着上面的例子，首先我们得根据三个自变量的范围进行采样，这边采样的方法一般都是蒙特卡洛采样以及一系列基于蒙特卡洛采样的变种，这个例子中我们采用了Sobolsequence，具体的采样原理在这就不说了，大家可以自行谷歌。为了方便讲解例子我们设置采样的样本数为4(N=4)，自变量数目为3(D=3)。我们按照上述网页的步骤。4.生成N*2D（即4行6列）的样本矩阵。这个就是我们Sobolsequence做的事情。这边我们生成的矩阵为：5、将矩阵的前D列设置为矩阵A，后D列设置为B列，在我们的例子中就是矩阵m的前3列设置为矩阵A，后3列设置为矩阵B。0.50.50.50.750.250.250.250.750.750.3750.3750.625A0.50.50.50.250.750.750.750.250.250.8750.3750.125B构造N*D的矩阵ABi(i=1,2,…,D)，即用矩阵B中的第i列替换矩阵A的第i列，以本体为例：241231Y=sin()+7(sin())+0.1sin()xxxx0.50.50.50.50.50.50.750.250.250.250.750.750.250.750.750.750.250.250.3750.3750.6250.8750.3750.125m10.50.50.50.250.250.250.750.750.750.8750.3750.625AB20.50.50.50.750.750.250.250.250.750.3750.3750.625AB30.50.50.50.750.250.750.250.750.250.3750.3750.125AB6.经过这三步我们构造了A、B、AB1、AB2、AB3这五个矩阵，这样我们就有(D+2)*N(即20)组x1、x2、x3输入数据，因此我们将有20组Y值。将上述的数据带入函，这里详细的计算过程就不描述了。根据输入我们得出对应的Y值矩阵。2.091363877681.110366058553.507651769321.31095036329AY2.091363877683.507651769321.110366058551.70665120003BY12.091363877680.675961634813.955626030581.71834424557ABY22.091363877683.934324819330.6836930085371.31095036329ABY32.091363877681.13166726983.499920395591.30537043023ABY7.根据一介影响指数公式：i(Y)()iXiXiSVarEXVarY总效应指数:i((Y))()TiXXiiSEVarXVarY其中:i11(Y)(B)((AB)(A))NiXiXijjjjVarEXfffN2i11((Y))((A)(AB))2NiXXiijjjEVarXffN(Y)(YY)ABVarVar这里的“加”不是普通意义上的相加而是构造了一个新的矩阵包括YA和YB即241231Y=sin()+7(sin())+0.1sin()xxxx2.091363877681.110366058553.507651769321.310950363292.091363877683.507651769321.110366058551.70665120003Y这样我们就可以通过上述的公式求解x1的一介影响指数。下面是求解过程。8.Y的均值等于mean(Y)=(2.09136387768+1.11036605855+3.50765176932+1.31095036329+2.09136387768+3.50765176932+1.11036605855+1.70665120003)/8=2.0545456218Var(Y)=0.8353325815421111()NXBjAjABjjVarYYYN本例中VarX1就等于(2.09136387768*(2.09136387768-2.09136387768)+3.50765176932*(0.67596163481-1.11036605855)+1.11036605855*(3.95562603058-3.50765176932)+1.70665120003*(1.71834424557-1.31095036329))/4=-0.0827611931562那么10.08276119315620.099075730.835332581542S12i11()2NXAjABjjEYYN本例中EX~1就等于((2.09136387768-2.09136387768)^2+(1.11036605855-0.67596163481)^2+(3.50765176932-3.95562603058)^2+(1.31095036329-1.71834424557)^2)/(2*4)=0.06941973968，那么10.069419739680.08310431220.835332581542TS至此x1的一介影响指数和全局影响指数都可以求出来了，然后x2、x3的一介影响指数与全局影响指数的求解过程和上面一样。当然在实际的求解过程中肯定不会用这么少的样本量的，但是求解过程是一样的。希望这篇文章会对做敏感性分析的小伙伴们有些帮助。