【MATLAB第111期】基于MATLAB的sobol全局敏感性分析方法二阶指数计算

一、简介

在MATLAB中计算Sobol二阶效应指数通常涉及到全局敏感性分析（Global Sensitivity Analysis, GSA），其中Sobol方法是一种流行的技术，用于评估模型输入参数的敏感性。Sobol二阶效应指数衡量的是两个参数之间的交互作用对模型输出的影响。

Sobol二阶效应指数的计算涉及到以下步骤：

1、生成Sobol序列，并在一阶的基础上，生成N=(2D+2)*npop个样本集，其中D为变量数， npop为采样的数量， N为总样本数 。

在一阶基础上，N=(D+2)*nPop样本， 包含A、AB和B矩阵。

二阶需要生成BA矩阵，用来评估二阶指数。

2、模型计算

可参考CSDN64期文章， 利用sobol函数进行抽样，得到的X值 ，通过bp组成的代理模型进行计算。

3、计算Sobol指数：使用Sobol序列和模型输出，计算每个参数的一阶、二阶效应指数和总效应指数， 其中，一阶和总效应指数较为好计算， 二阶效应指数可以参考python的Salib库进行研究 ，

Sobol二阶效应指数的计算公式如下：

其中i为1：D，j为i+1：D，

4、计算效果

二、部分源码

S2计算代码如下：

 Vjk = mean(BAj .* ABk - A .* B) / var(y);
  Sj = first_order(A, ABj, B);
  Sk = first_order(A, ABk, B);
  S2 = Vjk - Sj - Sk;

核心参考代码如下（需要自行二次编译）：

  % Normalize the model output
  Y = (Y - mean(Y)) / std(Y);
  % Separate output values
  A = Y(1:2*D+2:end);
  B = Y((end-1):-(2*D+1):1);
  AB = zeros(length(Y)/ (2*D+2), D);
    BA = zeros(length(Y)/ (2*D+2), D);
    for j = 1:D
      AB(:, j) = Y((j+1):2*D+2:end);
      BA(:, j) = Y((j+1+D):2*D+2:end);
    end
  end

  % Calculate second order indices if required

    for j = 1:D
      for k = j+1:D
        Si.S2(j, k) = second_order(A, AB(:, j), AB(:, k), BA(:, j), B);
end
end
function S = first_order(A, AB, B)
  % First order estimator following Saltelli et al. 2010 CPC, normalized by
  % sample variance
  y = [A, B];
  if range(y) == 0
    warning('Constant values encountered, indicating model evaluations (or subset of evaluations) produced identical values.');
    S = 0;
    return;
  end
  S = mean(B .* (AB - A)) / var(y);
end

function S = total_order(A, AB, B)
  % Total order estimator following Saltelli et al. 2010 CPC, normalized by
  % sample variance
  y = [A, B];
  if range(y) == 0
    warning('Constant values encountered, indicating model evaluations (or subset of evaluations) produced identical values.');
    S = 0;
    return;
  end
  S = 0.5 * mean((A - AB) .^ 2) / var(y);
end

function S = second_order(A, ABj, ABk, BAj, B)
  % Second order estimator f ollowing Saltelli 2002
  y = [A, B];
  if range(y) == 0
    warning('Constant values encountered, indicating model evaluations (or subset of evaluations) produced identical values.');
    S = 0;
    return;
  end
  Vjk = mean(BAj .* ABk - A .* B) / var(y);
  Sj = first_order(A, ABj, B);
  Sk = first_order(A, ABk, B);
  S = Vjk - Sj - Sk;
end