Sobol全局敏感性分析

参考敏感性分析—Sobol这篇文章，用python编写了一阶影响指数

import numpy
from sobol import sobol_seq
def fun(x1, x2, x3):
    Yvalue = numpy.sin(x1) + 7*(numpy.sin(x2))**2+0.1*x3**4*numpy.sin(x1)
    return Yvalue
def First_order(M, nchrom, lim_b):
    nchrom2 = nchrom*2
    seed = 2
    #result = sobol_seq.i4_sobol_generate(nchrom2, M, seed)
    # print result
    result = numpy.array([[0.5, 0.75, 0.25, 0.375], [0.5, 0.25, 0.75, 0.375], [0.5, 0.25, 0.75, 0.625], [
                         0.5, 0.25, 0.75, 0.875], [0.5, 0.75, 0.25, 0.375], [0.5, 0.75, 0.25, 0.125]])
    A = result[:nchrom]
    B = result[nchrom:]
    tmp_A = numpy.copy(A)
    AB = [0 for _ in range(nchrom)]
    for i in range(nchrom):
        tmp_A[i] = B[i]
        AB[i] = numpy.copy(tmp_A)
    YA = fun(A[0], A[1], A[2])
    YB = fun(B[0], B[1], B[2])
    YAB = [0 for _ in range(nchrom)]
    for i in range(nchrom):
        YAB[i] = fun(AB[i][0], AB[i][1], AB[i][2])
    VarAB = [0 for _ in range(nchrom)]
    for j in range(nchrom):
        sigmaAB = 0
        for i in range(M):
            sigmaAB = sigmaAB + YB[i]*(YAB[j][i]-YA[i])
        VarAB[j] = sigmaAB/M
    Y = YA, YB
    Y = numpy.reshape(Y, (8, 1))
    imp_index = [0 for _ in range(nchrom)]
    for i in range(nchrom):
        imp_index[i] = VarAB[i]/numpy.var(Y)
    return imp_index
nchrom = 3
M = 4
lim_b = [[0, 1], [1, 2], [3, 4]]
imp_index = First_order(M, nchrom, lim_b)
print imp_index

写完之后发现有第三方包。