數據污染情形下的全局靈敏度分析

謝 恩，馬義中，劉麗君，林成龍

(南京理工大學 經濟管理學院，江蘇 南京 210094)

0 引言

全局靈敏度分析通常假設試驗數據服從正態分布，但數據污染使得實際試驗數據偏離假設,伴隨著大數據時代海量數據的出現，數據被污染的風險加大[1]。同樣，基于污染數據的仿真試驗可能會歪曲事物本來的面目，降低說服力，甚至可能得出謬誤的結論。因此，如何改進統計方法，使得在面對數據污染時仍能夠得到穩健的結果，成為近年來研究的熱點問題[2-3]。樣本中存在異常值被認為是數據污染的一種，傳統的異常值檢驗技術[4]通過識別異常值并直接舍棄的處理方法，往往導致表示過程質量特性的重要信息丟失，影響后續工作[5]。RIPLEY[6]指出在多元或高維數據情形下難以識別數據中的異常值，需要采用穩健統計方法來抵抗污染值的干擾。

基于穩健模型和穩健統計技術建立模型并估計模型參數，能更好地描述數據集中大多數數據的統計特性，能抵御異常值對仿真試驗的影響[7]。SERFLING[8]指出樣本數據存在污染的情形時，用中位數和中位數絕對偏差(Median Absolute Deviation, MAD)代替均值和方差估計樣本位置和尺度參數可以獲得穩健性較高的結果。韓云霞等[9]研究了基于穩健統計量(Hodges-Lehmann, HL)和Shamos估計數據污染情形下最優參數置信區間，結果表明，基于穩健統計量的方法估計精度更高，同時能夠抵抗污染值的干擾。PARK等[10]研究了數據污染情形下不同的穩健統計量估計位置參數和尺度參數，通過分析蒙特卡洛仿真試驗結果，給出幾種穩健統計量在數據污染情形下的有效性，其中分別用HL和Shamos估計位置參數和尺度參數進行穩健設計的最優解優于其他穩健統計量。劉麗君等[11-12]用中位數估計位置參數，用MAD估計尺度參數，改進了序貫分支試驗的顯著性檢驗程序，在不增加額外試驗次數的前提下，解決了不同類型數據污染情形下的因子篩選問題。

靈敏度分析方法主要分為3類[13-14]：①旨在區分模型輸入因子是否對模型輸出有影響的因子篩選方法(screening method)；②局部靈敏度分析方法(local sensitivity analysis)，考慮圍繞特定設計點的微小擾動對模型輸出不確定性影響；③全局靈敏度分析方法(global sensitivity analysis)，考慮整個輸入空間內的因子變化對模型輸出不確定性的影響，旨在對模型所有輸入因子基于其對模型輸出不確定效應的貢獻進行重要度排序。Sobol’指數法作為基于方差的全局靈敏度分析方法[15-16]，獨立于模型輸入和輸出，容易解釋和實現，適用于模型輸入因子排序和重要變量篩選，在工業工程、環境和大氣工程及化工工程等領域被廣泛研究和使用[17-19]。但在面對復雜問題時，Sobol’指數通過大量的模型估計也很難得到收斂的合理解，無法識別因子重要度。隨著計算機技術的發展和進步，使得幾乎所有工程領域復雜系統的行為都能用數值模型來模擬和預測[20]，因此，基于代理模型的全局靈敏度分析方法獲得了高度關注[20-22]。另一方面,基于方差的Sobol’指數蒙特卡洛仿真計算技術也得到了改進[23-25]，雙循環重排序方法(Double Loop Reordering approach, DLR)作為Sobol’蒙特卡洛仿真技術改進方法的一種，使得仿真計算效率得到很大的提高[26-29]。

然而，在數據存在污染的情形下，不能準確反映樣本的統計特性，導致基于方差的靈敏度分析方法無法識別因子重要度或對因子重要度進行排序，無法正確度量模型輸出不確定性的來源，從而使得模型更加復雜、計算成本更高。針對仿真試驗中數據污染(本文指數據中存在異常值)的問題，引入穩健統計量改進Sobol’法中的DLR蒙特卡洛仿真方法，提出了穩健雙循環重排序方法(Robust Double Loop Reordering approach, RDLR)，并通過仿真試驗驗證了所提方法的有效性和穩健性。

1 穩健統計量及其性質

1.1 穩健統計量及其相對效率

1.2 兩組穩健統計量的漸近分布特性

基于穩健統計量HL-Shamos組合，可構造如下統計量[30]：

(1)

基于穩健統計量Med-MAD組合，可構造如下統計量[33]：

(2)

為了驗證兩組穩健統計量的分布情況，在正態分布假設條件下，分別在樣本量(n)為10和100時進行1 000次隨機抽樣，基于樣本點繪制統計量TMM，THS的正態分布QQ圖，以及樣本量為10的經驗累計分布圖和概率密度函數圖如圖1所示。

對比圖1a和圖1b可知，在樣本量較大時，兩組統計量更趨向于正態分布；在樣本量較小時，基于統計量THS的分布更趨向于正態分布。如圖1c所示為基于統計量THS、TMM以及均值—方差的經驗累計分布圖，如圖1d所示為基于統計量THS、TMM以及均值—方差的密度函數，分析圖1c和圖1d可得，統計量THS、TMM的分布均近似服從正態分布，且基于統計量THS的分布更接近正態分布。因此，在樣本量較小時建議優先選擇THS統計量?；诜€健統計量的漸近正態分布特性，本文采用穩健統計量代替均值和標準差估計樣本的位置參數和尺度參數是合理的。

2 基于穩健統計量的全局靈敏度分析

2.1 穩健統計量替代均值和方差的全局靈敏度分析

假設系統某一質量特性Y和m個可控輸入變量x=(x1,…,xm)之間的關系為Y=f(x)，則模型輸出的全方差公式為：

V[Y]=Vxi[Ex-i[Y|xi]]+Exi[Vx-i[Y|xi]]。

其中：E[·]表示期望，V[·]表示方差，x-i表示除xi之外的輸入因子，Vxi[Ex-i[Y|xi]]用來定量地度量輸入因子xi對模型輸出的影響。因此，因子xi基于方差的一階全局靈敏度指數計算如下：

Si=Vxi[Ex-i[Y|xi]]/V[Y]。

當模型輸出中含有異常值時，考慮采用穩健統計量HL-Shamos或者Med-MAD替代均值—方差計算全局靈敏度指數，使得基于穩健統計量的全局靈敏度分析方法在面對數據污染時仍能準確識別模型輸入因子的重要度。

其中：a=1.158 75;a1=0.414 253 297;a2=0.442 396 799;當n≤100時，a3=a4=0，當n>100時,a3=2.822，a4=12.238。

其中：b=2.702 7;b1=-0.762 13;b2=-0.864 13;當n≤100時，b3=b4=0,當n>100，且n為奇數時，b3=0.299 6，b4=-149.357，當n為偶數時,b3=-2.417，b4=-153.01。

基于穩健的位置參數—尺度參數為HL-Shamos組合的全局靈敏度指數定義為Si1,i2,…,ik(HS)，

Si1,i2,…,ik(HS)=

(3)

基于穩健的位置參數—尺度參數為Med-MAD組合的全局靈敏度指數定義為Si1,i2,…,ik(MM)，

Si1,i2,…,ik(MM)=

(4)

式(3)和式(4)中{i1,i2,…,ik}是{1,2,…,m}的一個子集，且1≤i1<…

2.2 基于穩健統計量的RDLR蒙特卡洛仿真計算方法

假設函數f在積分空間Im=[0,1]m上平方可積，函數f通過Sobol’分解可分解為：

fi,j(xi,xj)+…+f1,2,…,m(x1,x2,…,xm)。

當文獻[15]中的條件滿足時，上式的分解是唯一的，對其兩邊平方并積分，可得：

假定兩互補子集y和z構成輸入變量x=(y,z)，令子集y=(xi1,xi2,…,xik)，子集y的方差

Sobol’法的總方差為

其中{i1,i2,…,ik}是{1,2,…,m}的一個子集，即1≤i1<…

在仿真試驗中，假定x和x′是樣本空間中兩個N×m維相互獨立的抽樣數組，令x=(y,z)，x′=(y′,z′)，對應于輸入因子子集y=(xi1,…,xik)的方差為：

針對樣本均值和方差對異常值敏感的問題，提出穩健雙循環重排序方法(RDLR)，采用穩健統計量替代均值—方差改進傳統的DLR方法。用穩健的位置統計量估計每個分區中Nm個模型輸出的位置參數，即對第k(1≤k≤M)個分區分別計算HL統計量和中位數Med統計量來估計每個分區中模型條件輸出的位置參數，形式如下：

p,q=1,…,Nm,k=1,…,M；

由上式可得M個穩健的條件輸出位置統計量，然后用穩健的尺度參數統計量Shamos和MAD估計M個位置參數統計量的穩健尺度參數：

(5)

(6)

模型非條件輸出的穩健尺度參數為：

f(yq,zq)|))2/A(N),p,q=1,…,N；

(7)

(8)

因此，由式(5)～式(8)可得，輸入變量y=xi的基于穩健統計量改進的RDLR方法的一階全局靈敏度指數，即一階Sobol’指數的計算公式為：

Sy(HS)=Dy(HS)/DHS；

(9)

Sy(MM)=Dy(MM)/DMM。

(10)

2.3 RDLR方法的實現步驟

基于穩健統計量的RDLR方法具有較高的穩健性，能夠抵御異常值的影響，其仿真算法實現步驟如下：

步驟1算法初始化(測試函數選擇)，隨機生成N個m維的樣本點xj，計算對應的輸出f(xj)。

步驟2令y=xi(1≤i≤m)獲得排序后的樣本點x(j),j=1,2,…,N，將排序后的輸入變量對應的輸出f(xj)分成M(M

步驟4計算所有樣本點xj對應的非條件輸出的穩健尺度參數DHS和DMM。

步驟5計算基于穩健統計量的輸入因子靈敏度指數Dy(HS)/DHS和Dy(MM)/DMM。

步驟6重復步驟2～步驟5，令y取每一個輸入變量，計算所有輸入變量的穩健靈敏度指數。

傳統的DLR和改進的RDLR方法能計算單個輸入變量的靈敏度指數，對變量子組無效。

3 算例分析

為了更好地說明所提方法的有效性，選用3個測試函數對提出的RDLR方法進行驗證。Ishigami函數和Sobol G—函數都是全局靈敏度分析研究中廣泛采用的測試函數。Ishigami函數具有輸入變量間相互關系復雜，能夠代表大多數模型輸入變量間關系類型的特點[34]；Sobol G—函數輸入變量數目可變，并且可以通過調節函數表達式中非負常數，進而達到改變模型輸入變量的靈敏度指數的目的，具有較高的靈活性[35]?；诜讲畹姆椒o法識別非線性偏態分布輸出函數的輸入變量的重要度[36]，為了驗證本文基于穩健統計量的方法針對非正態問題的有效性，選擇該二維輸入的非線性偏態分布輸出函數。

考慮數據污染對基于方差的全局靈敏度分析方法的影響，人為增加一個初始輸出的污染函數，記為fcon(y,λ)，污染后的模型輸出為：fcon(y,λ)=y+m·s·I(y),數據污染的相關函數λ=(p,m,s)與3個參數有關，其中p表示樣本數據中變異數據的概率，本文取值分別為0(不含異常值)、0.05和0.1，此處給定污染數據占比均小于各統計量的失效點；m表示變異量，取值分別為100和300；s表示變異符號，取值為1和-1；其中I(y)為示性函數，即I(y≥0)=1，I(y<0)=-1?；诟鹘o定變異函數污染參數的取值，可得9種污染參數組合，分別命名為Ⅰ、Ⅱ、Ⅲ、Ⅳ、Ⅴ、Ⅵ、Ⅶ、Ⅷ和Ⅸ，第Ⅰ種變異概率為0，污染參數組合為(0,-,-)，模型輸出無污染,不同污染情形的參數組合如表1所示。

表1 不同數據污染情形的參數組合

3.1 Ishigami函數

Ishigami函數[34]包含3個輸入因子，表達式為：

f(x1,x2,x3)=sinx1+a·sin2x2+

其總方差D和條件方差D1,D2,D3的理論計算值如下所示：

令a=7，b=0.1，Ishigami函數基于方差的輸入因子x1,x2,x3的全局靈敏度指數解析解分別為：S1=0.313 8,S2=0.442 4,S3=0，因此S2>S1>S3,即因子x2重要度為第一，因子x1重要度為第二，因子x3重要度為第三。

如圖2所示為Ishigami函數1 000次仿真試驗輸出在9種不同的污染參數組合情形的均值、HL、中位數、標準差、Shamos和MAD等6個統計量的條形圖。

如圖2a所示為不同的位置參數統計量的條形圖，對比分析可知樣本均值對異常值非常敏感，而HL和中位數幾乎不受異常值影響；圖2b所示為不同尺度參數統計量的條形圖，其中標準差受到異常值的影響最大，穩健尺度統計量Shamos和MAD受異常值影響較小。因此，面對數據污染問題，基于穩健統計量的仿真試驗結果不易受異常值干擾。

如表2所示為Ishigami函數輸出在9組污染情形下的不同統計量的方差，受異常值影響，9組輸出的均值和標準差的方差遠遠大于穩健統計量的方差，即均值和標準差對異常值更加敏感。

表2 Ishigami函數輸出在9種污染參數情形下不同統計量的方差

為驗證數據污染情形下基于穩健統計量的全局靈敏度分析方法的有效性和穩健性，本試驗分別在9種不同污染參數組合情形下各進行1 000次獨立重復試驗，得到因子靈敏度指數。在每種污染參數組合條件下、每次獨立試驗中，在模型輸入空間中隨機抽取1 000個樣本點，計算得到1 000個輸出，執行一次靈敏度指數計算，得到一組靈敏度指數。Ishigami函數輸入因子靈敏度指數基于方差、HL/Shamos和Med/MAD等方法在不同污染情形下，試驗所得靈敏度指數的均值統計結果如表3所示。

由表3可知，基于方差的全局靈敏度分析方法，僅在第I種污染參數組合，即輸出無污染的情形下能識別因子的重要度；在輸出受到污染的情形下，因子x1的靈敏度指數S1和因子x2的靈敏度指數S2幾乎相等，無法識別輸入因子重要度?；贖L/Shamos和Med/MAD的全局靈敏度指數在9種污染參數組合情形下，都能識別輸入因子的重要度，可得輸入因子靈敏度指數關系為S2>S1>S3，與解析解的結果一致；對比基于兩組不同穩健統計量的因子靈敏度指數仿真試驗結果，基于HL/Shamos方法的性能更優，各因子靈敏度指數更接近解析解。

表3 基于不同統計量的Ishigami函數在不同污染情形下因子靈敏度指數

如圖3所示為9種污染情形下，Ishigami函數基于不同統計量各進行1 000次獨立重復試驗的因子靈敏度指數箱線圖。由圖3a可知，基于方差的方法在數據污染的情形下，不能有效識別因子重要度，基于方差方法對異常值敏感，仿真試驗結果受到異常值影響。圖3b和圖3c分別為基于穩健統計量HL/Shamos和Med/MAD方法的靈敏度指數箱線圖，結果表明：在9種污染情形下基于穩健統計量的全局靈敏度分析方法都能正確識別因子重要度，因子靈敏度指數關系為S2>S1>S3，與解析解的結果一致。因此，基于穩健統計量的靈敏度指數計算方法能夠實現數據污染情形下因子重要度排序，相比基于方差的方法更穩健。比較圖3b與圖3c中第Ⅰ種污染參數組合情況可知在輸出無污染時，基于穩健統計量的全局靈敏度分析方法同樣能夠正確識別輸入因子的重要度。

樣本量較小時，由于統計量THS相對于統計量TMM的分布更趨向正態分布，且統計量HL/Shamos相對統計量Med/MAD的效率更高，因此基于統計量HL/Shamos的全局靈敏度分析方法性能更優。

3.2 非線性偏態分布輸出函數

輸出高度偏態分布的函數[36]：

y=x1/x2,x1～χ2(10),x2～χ2(13.978)。

輸入變量x1對函數輸出的貢獻比x2對函數輸出的貢獻大，然而基于方差的全局靈敏度分析方法無法識別出兩因子的重要度。

為了驗證輸出偏態分布且受到污染情形下基于方差、HL/Shamos和Med/MAD等方法的有效性和穩健性，采用與3.1節相同的試驗程序，計算9種污染情形下基于不同統計量的輸入因子靈敏度指數。在顯著水平為0.05時，對9種污染條件下基于方差方法的兩因子1 000次仿真試驗的靈敏度指數進行檢驗。結果顯示，除第I種污染參數組合外，其它污染參數組合下的P值都大于顯著水平，因此兩因子1 000次仿真試驗的靈敏度指數無顯著差異，不能對輸入因子進行重要度排序，t檢驗P值如表4所示。

表4 不同污染情形下基于方差方法的因子靈敏度指數檢驗

如表5所示為輸出偏態分布的函數在9種污染情形下基于不同統計量的全局靈敏度分析方法，各進行1 000次獨立重復試驗，所得因子的靈敏度指數均值統計結果。由表5可知，在函數輸出偏態分布或被污染情形時，基于方差全局靈敏度分析方法無法識別因子重要度，即因子x1的靈敏度指數S1和因子x2的靈敏度指數S2的關系為S1≈S2。對比基于穩健統計量HL/Shamos和Med/MAD方法的靈敏度指數可知，因子x1的靈敏度指數S1顯著大于因子x2的靈敏度指數S2，因此基于穩健統計量的靈敏度分析方法不受數據污染或偏態分布的影響，能正確識別因子重要度，具有較高的穩健性。

表5 不同污染情形下基于不同統計量的因子靈敏度指數

如圖4所示為輸出偏態分布的函數在9種不同污染情形下，基于不同的統計量的全局靈敏度分析方法，各進行1 000次獨立重復試驗，所得因子靈敏度指數的箱線圖。

由圖4a可知，在函數輸出呈偏態分布條件下，不論輸出是否存在污染的情形，基于方差的方法均無法識別因子的重要度，仿真試驗結果無法對因子重要度排序；圖4b和圖4c是基于穩健統計量的靈敏度指數箱線圖，結果顯示函數輸出在9種污染情形下，輸入因子x1的全局靈敏度指數S1顯著大于輸入因子x2的全局靈敏度指數S2，因此基于穩健統計量的全局靈敏度分析方法在數據偏態分布情形下，不論數據是否被污染都能準確識別輸入因子的重要度；基于穩健統計量的全局靈敏度分析具有較高的穩健性和更廣泛的適用性。

3.3 Sobol G-函數

在較高維輸入和較高污染數據比例情形下，為了驗證本文基于穩健統計量的全局靈敏度分析方法的有效性，選用全局靈敏度分析中常用的測試函數Sobol G—函數[35]，其形式為：

其中：輸入因子xi(i=1,2,…,d)在[-1,1]d上均勻分布，ai是非負常數。Sobol G—函數可以通過控制非負常數來調節每一個輸入對輸入方差的貢獻。當非負常數ai小時，其對應的輸入因子對輸出的方差貢獻就更大;當非負常數ai大時，其對應的輸入因子對輸出的方差貢獻就小。本文設置輸入因子數目為6，每個輸入因子對應的非負常數是(a1,a2,a3,a4,a5,a6)=(15,3,2,1,0.1,0)，模型輸出的污染參數組合如表6所示。

表6 Sobol G—函數的輸出污染參數組合

考慮穩健統計量HL和Shamos的極限BP點為29.3%，在此驗證數據中污染數據占比較高的情形下，基于穩健統計量的全局靈敏度方法的有效性，設置污染參數最大比例為25%。如表7所示為不同污染情形下6維Sobol G—函數基于不同統計量的1 000次仿真試驗所得輸入變量靈敏度指數均值統計結果。

表7 不同污染情形下Sobol G—函數基于不同統計量的因子靈敏度指數

由表7可知，基于方差的方法僅在數據不存在污染的情形下能夠正確地對輸入因子重要度進行排序，再次驗證了基于方差的全局靈敏度分析方法不適用于數據污染的情形?；诜€健統計量的方法既能在數據無污染情形下識別因子重要度，又能在數據存在污染的情形下正確地識別因子的重要度。當數據中污染值的比例達到25%時，基于穩健統計量Med/MAD的方法優于基于穩健統計量HL/Shamos的方法，可能的原因是穩健統計量Med/Mad的BP點較高，而25%的污染值的占比使得穩健統計量HL/Shamos的性能較低。

4 結束語

針對數據污染情形下基于方差的全局靈敏度分析方法不能正確識別模型輸入因子重要度的問題。本文提出一種新的Sobol’指數蒙特卡洛仿真計算方法RDLR法。通過蒙特卡洛仿真計算兩個測試函數在不同的污染情形下基于不同方法的因子靈敏度指數，對比傳統的DLR方法和本文所提方法的穩健性和有效性。分析試驗結果可得出以下結論：

(1)本文所提方法具有良好的穩健性，在模型輸出存在異常值時，能正確識別因子重要度，有效抵御異常值的影響；

(2)RDLR法具有廣泛的適用性，在模型輸出理想無污染或高度非正態分布的情形下，能正確識別因子重要度；

(3)RDLR方法不必增加仿真試驗次數，保證仿真試驗的經濟性；

(4)在模型輸入因子數目較小時，基于穩健統計量HL/Shamos的方法性能更優，當數據中污染數據占比較高的情形下，基于穩健統計量HL/Shamos的方法性能更優。

本文主要貢獻是通過引入穩健統計量改進了Sobol’蒙特卡洛仿真計算方法，進而提高了其抗異常值特性，拓展了全局靈敏度分析方法的適用范圍。需要指出的是，本文假設模型為單個輸出并且輸入因子相互獨立，當模型有多個輸出或輸入因子具有相關關系時，穩健的全局靈敏度分析方法將是一個值得關注的方向。