敏感性分析

您所在的位置:网站首页 灵敏度分析举例 敏感性分析

敏感性分析

2023-08-08 23:27| 来源: 网络整理| 查看: 265

最近在搞敏感性分析方法中的Sobol方法,看了一些国内的论文,发现一个通病,就是公式一挂就可以得出结果了,对于我这种数学渣渣真的很不适合。直到我在维基百科上面找到了这种方法的详细解释https://en.wikipedia.org/wiki/Variance-based_sensitivity_analysis,包括分析的步骤都写的很清楚,当然还有蛮多东西自己还不是太了解(越来越发觉,数学真的太重要了,师弟师妹们大学别光玩了,至少工科数学还是要学好的。)。本篇博客就不对概念性的东西再做阐述,主要自己目前也是半吊子,就不在这误导大家了。 分析的步骤刚才的网页已经有详细的描述,今天我们就根据网页上的步骤用一个例子来走一遍。 假设现在有一个函数 这里写图片描述,从公式中可以看到有x1、x2、x3三个自变量对应变量Y有影响。然后一般会给这三个参数一个取值范围,这里假设三个自变量的取值范围都设为[0,1]。敏感性分析的目的就是求取这三个参数对于Y值得贡献。当然我们这边可能有人一下子就可以分析出那个参数对于Y值影响最大,但是在解决实际问题时,这个函数一般都是未知,我们只能将其视作一个黑盒子,只有输入和输出,这时我们对其进行敏感性分析就很有必要了。经过敏感性分析我们就能找出对结果影响较大的参数。这样对于调整结果是很有帮助的。 接着上面的例子,首先我们得根据三个自变量的范围进行采样,这边采样的方法一般都是蒙特卡洛采样以及一系列基于蒙特卡洛采样的变种,这个例子中我们采用了 Sobol sequence,具体的采样原理在这就不说了,大家可以自行谷歌。为了方便讲解例子我们设置采样的样本数为4(N = 4),自变量数目为3(D = 3)。我们按照上述网页的步骤。 1、 生成N * 2D(即4行6列)的样本矩阵。这个就是我们Sobol sequence做的事情。这边我们生成的矩阵为:

这里写图片描述 2、 将矩阵的前D列设置为矩阵A,后D列设置为B列,在我们的例子中就是矩阵m的前3列设置为矩阵A,后3列设置为矩阵B。 这里写图片描述 构造N*D的矩阵ABi(i = 1,2,…,D),即用矩阵B中的第i列替换矩阵A的第i列,以本体为例: 这里写图片描述 经过这三步我们构造了A、B、AB1、AB2、AB3这五个矩阵,这样我们就有(D + 2) * N (即20)组x1、x2、x3输入数据,因此我们将有20组Y值。将上述的数据带入函数 这里写图片描述,这里详细的计算过程就不描述了。根据输入我们得出对应的Y值矩阵。 这里写图片描述这里写图片描述 4、 根据一介影响指数公式: 这里写图片描述 总效应指数:这里写图片描述 其中这里写图片描述 这里写图片描述 这里写图片描述 这里的“加”不是普通意义上的相加而是构造了一个新的矩阵包括YA 和 YB 即这里写图片描述 这样我们就可以通过上述的公式求解x_1的一介影响指数。下面是求解过程。 Y 的均值等于 mean(Y)=(2.09136387768+1.11036605855+3.50765176932+ 1.31095036329+2.09136387768+3.50765176932+ 1.11036605855+1.70665120003)/8 = 2.0545456218 Var(Y) = 0.835332581542 这里写图片描述 本例中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 那么 这里写图片描述 这里写图片描述 本例中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,那么这里写图片描述 至此x1的一介影响指数和全局影响指数都可以求出来了,然后x2、x3的一介影响指数与全局影响指数的求解过程和上面一样。 当然在实际的求解过程中肯定不会用这么少的样本量的,但是求解过程是一样的。希望这篇文章会对做敏感性分析的小伙伴们有些帮助。



【本文地址】


今日新闻


推荐新闻


    CopyRight 2018-2019 办公设备维修网 版权所有 豫ICP备15022753号-3