一種求解相關隨機變量聯合累積分布函數的方法

楊 永，李海濱

(內蒙古工業大學 理學院，內蒙古自治區 呼和浩特 010051)

聯合概率累積分布函數是概率論中重要的概念，在工程結構可靠性分析中有著廣泛應用.實際應用中，根據已有的一維函數模型，如正態分布、對數正態分布、極值Ι型分布、指數分布、威布爾分布等，能夠容易得到邊緣累積分布函數，對具有獨立性的數據，聯合累積分布函數可以通過邊緣累積分布函數的直接相乘得到.但在工程中，所測得的數據往往存在相關性，如巖土工程抗剪強度參數中的黏聚力和內摩擦角、基樁荷載—位移雙曲線等.為了能更好的進行可靠性分析，必須考慮數據間的相關性.

目前，為了計算方便，針對具有相關性的二維數據往往采用二維分布模型的方法，如二維正態分布[1]、二維對數正態分布[2]、二維指數分布[3]等.雖然這些二維分布模型能夠考慮到數據間的相關性，但是必須滿足邊緣分布具有相同的分布類型，如二維正態分布要求其邊緣分布均滿足正態分布.但在實際工程中，數據往往不是同類型分布，因此二維分布模型的適用范圍受到了極大的限制.對于這些非同類型分布，Nataf于1962年提出了Nataf分布模型[4]，有效地解決了一些相關非正態變量的聯合概率分布問題，但是基于Nataf模型構造的二維聯合累積分布函數依舊隱含著變量間的相關模型是Gaussian結構，因此還是很難解決工程問題中復雜邊緣累積分布函數的問題.另外一種常用的構造方法是基于Pearson以及Spearman相關系數構造出聯合累積分布函數，但是其只能描述數據間的線性相關關系，很難處理一些具有明顯非線性特征的問題.針對現有構造聯合累積分布函數存在的問題，尋找一種適用范圍廣、計算簡便的新方法是一個亟需解決的問題.

1959年，Sklar指出，任意一個多維聯合累積分布函數都可以分解為對應的邊緣累積分布函數和一個Copula函數[4].Copula理論的提出為聯合累積分布函數的構造提供了一個全新的途徑，由于分開構造邊緣累積分布函數和Copula函數，因此可以構造出任意邊緣累積分布函數和任意相關結構的聯合累積分布函數，而且對不同相關結構的數據，有多種Copula函數可以選擇，如Gaussian Copula、t-Copula、Frank Copula等.由于Copula函數的優良性質，其在金融、水文、巖土等多個領域得到了廣泛應用[5].而阿基米德Copula函數是所有Copula函數中最常用的一類函數，具有形式簡單、對稱性、可結合性等優點，而且只要找到“生成元”就可以構造出這一類函數.

目前已經有許多構造阿基米德Copula函數的方法，如Markov方法[6]、Laplace變換[7-8]等，但這些方法均有一定的局限性[9]，因此一些學者對阿基米德Copula生成元的復合構造進行了研究[10].在此基礎上，文中提出了一種多參數Copula構造形式，并利用經驗分布函數，通過采用正交距離回歸算法進行曲線擬合的方法，對由此函數構造的聯合累積分布函數進行擬合，從而實現了對聯合累積分布函數的求解.

首先介紹了最常用的阿基米德Copula函數及其構造具有相關性的二維隨機變量聯合累積分布函數的方法，然后提出一種多參數Copula函數的構造方法以及對其構造的聯合累積分布函數進行擬合和參數估計，最后利用算例對其有效性做出驗證.

1 基于Copula理論的聯合概率分布函數構造方法

Copula函數是對多元隨機變量相關結構的一種刻畫，可以理解成多元隨機變量的聯合累積分布函數和邊緣累積分布函數之間聯系的紐帶[9].

1.1 Copula函數的定義

Copula理論最早由Sklar于1959年提出.Sklar指出，任意一個多維聯合累積分布函數都可以分解為相應的邊緣累積分布函數和一個Copula函數，該Copula函數確定了變量間的相關性，包括相關系數的大小和相關結構的類型.對于n維情形,根據Sklar定理，可將變量x1,x2,…,xn的聯合累積分布函數F(x1,x2,…,xn)表示為：

F(x1,x2,…,xn)=C(F1(x1),F2(x2),…,Fn(xn),θ)=C(u1,u2,…un;θ),

(1)

式中，ui=Fi(xi)為變量xi的邊緣累積分布函數，同下文i=1，2，…,n；C(u1,u2,…un;θ)為Copula函數；θ為Copula函數的相關參數.

常用阿基米德Copula函數的具體形式為:

(2)

(3)

(4)

1.2 Copula函數的相關參數θ

Copula函數中的相關參數θ表征了變量間的相關性的大小，由于描述的是變量總體間的相關性，因此，需要借助極大似然理論來求出相關參數.用于對Copula函數進行參數估計的方法主要有最大似然估計法(ML估計)，分布估計(IMF估計)和半參數估計(CML估計)等[11-12].將采用半參數估計法也叫偽最大似然估計法進行參數估計，該方法通過采用樣本數據各個變量的經驗分布函數來分別取代其邊緣累積分布函數，以此得出Copula函數的相關參數θ.以具有邊緣分布F1(x1)和F2(x2)的二維聯合累積分布函數為例：

由(1)得樣本(x1i,x2i)，i=1,2,…，n的似然函數為:

(5)

其對數似然函數為:

(6)

對于(6)中的邊緣累積分布函數值，采用經驗分布函數值Femp代替，然后根據最大似然估計方法，求解出對數似然函數的最大值點，即可估計出相關參數θ.

(7)

1.3 Copula函數模型的擬合優度檢驗

判定選用的函數模型是否能夠較好的描述變量之間的關系，需要對Copula函數進行擬合優度檢驗，理論上，統計學常用的變量分布假設檢驗方法均適用于Copula函數檢驗，如皮爾遜擬合優度χ2檢驗、Kolmigrov檢驗、正態W檢驗等[13].

文中采用Kolmigrov-Smimov(K-S)來檢驗Copula聯合累積分布函數模型，用離差平方和準則(OLS)、赤池信息準則(AIC)對Copula聯合累積分布函數模型進行優選.

以二維為例，K-S檢驗統計量D、離差平方和準則(OLS)、赤池信息準則(AIC)定義如下：

1) K-S檢驗統計量D

(8)

其中，Ci是樣本數據xi=(x1i,x2i)的理論聯合累積分布函數值；mi是二維數據樣本數據中滿足條件x1≤x1i,x2≤x2i的個數；n為樣本數據量.

2) 離差平方和準則(OLS)

(9)

其中，

Femp(x1i,x2i)=p(x1≤x1i,x2≤x2i)，i=[1,n].

Femp(x1i,x2i)為經驗分布函數值，C(u1i,u2i)為理論聯合累積分布函數值.

3) 赤池信息準則(AIC)

(10)

AIC=nln(MSE)+2k,

其中，Femp(x1i,x2i)為經驗分布函數值，C(u1i,u2i)為理論聯合累積分布函數值，k為模型參數的個數.

2 一種乘積生成元構造的多參數Copula

不同的Copula函數用于描述具有不同相關特性的相關變量，在Copula函數中，阿基米德Copula函數是現在最常用也是最重要的一類Copula，由Genest和Mackey于1986年所提出，其函數表達式為：

C(u1,u2,…,un;θ)=φ-1(φ(u1)+φ(u2)+…+φ(un)),

(11)

式中φ(·)是阿基米德Copula的生成元，為滿足固定邊值的單調遞減凸函數.

根據生成元不同，阿基米德函數可以分為多種不同的形式，本文選取其中最常用的Gumbel Copula函數、Clayton Copula函數、Frank Copula函數等3種函數，其生成元和函數具體形式如表1.

表1 阿基米德Copula函數及其生成元

2.1 一種多參數Copula函數的構造

由上文可知，在阿基米德Copula函數中只有一個相關參數θ.由相關理論可知，生成元的乘積仍為生成元.以二維函數為例，將上述表1中的3種生成元兩兩相乘得到3種新的乘積生成元列于表2.在表2中給出了各乘積生成元所對應的隱式Copula函數，其中θ1、θ2為待定的相關參數.與已有的阿基米德Copula函數相比，新構造的Copula函數相關參數增加到2個，但由于選用了雙參數表示樣本的整體相關性，所以函數為不容易顯性化的隱函數，因此需要一種求解相關參數的新方法.

表2 乘積生成元及其隱式Copula函數

2.2 相關參數θ1、θ2的確定方法

曲線擬合是采用連續曲線去近似地刻畫一些離散點的函數關系，是處理數據常用的方法之一.主要包括三方面的內容，一是確定需要擬合的數據樣本點，二是選擇合適的函數模型，三是選取算法確定函數模型中的未知參數.在2.1節中，已經給出了邊緣累積分布函數與聯合累積分布函數之間含有未知參數的函數模型，因此，只要確定擬合數據以及算法就可以確定相關參數θ1、θ2.

2.2.1 擬合數據

在統計學中，經驗分布函數是對已有樣本累積分布函數的一種估計，根據Glivenko-Cantelli定理，隨著樣本數的增加，經驗分布逐步收斂于真實累積分布值，因此采用邊緣分布的經驗分布函數作為邊緣累積分布函數.而對各個樣本點處的聯合累積分布函數值，同樣可由聯合經驗分布函數代替.以前者為自變量，后者為因變量，即可構成擬合樣本.以n組二維數據樣本(x1,x2)為例：

邊緣累積分布函數

(12)

其中，

{x11,x12,…，x1n}*表示集合{x11,x12,…，x1n}中不大于x1i的個數,

{x21,x22,…，x2n}*表示集合{x21,x22,…，x2n}中不大于x2i的個數.

聯合累積分布函數：

(13)

其中，{(x11,x21),(x12,x22),…，(x1n,x2n)}*表示集合{(x11,x21),(x12,x22),…，(x1n,x2n)}中x1≤x1i且x2≤x2i數據點的個數.

2.2.2 正交距離曲線擬合

由于采用了經驗分布函數，為考慮邊緣經驗分布函數代邊緣累積分布函數所帶來的誤差，即自變量誤差，所以文中算法選取正交最小二乘法.正交最小二乘算法又稱正交距離曲線擬合，其與普通最小二乘的區別是考慮了函數擬合過程中自變量的誤差，以正交距離的殘差平方和極小為準則進行曲線擬合，使擬合結果從整體上達到擬合最佳[14].

正交距離回歸(ODR)算法[15]通過在迭代過程調整擬合參數使變量殘差的平方和最小.ODR中的殘差不是觀察值與變量的預測值之間的差異，而是從數據到擬合曲線的正交距離.擬合準則如下：

(14)

其中，f(xi+δxi,yi+δyi,β)=0i=1,……，n,wxi和wyi為自變量與因變量的權重系數，δxi、δyi為自變量xi和因變量yi的殘差，β是擬合參數.

3 數值算例

算例1

隨機生成100組數據點作為樣本.其中，隨機變量x1為正態分布，均值μ1=8，標準差σ1=2，隨機變量x2也為正態分布，均值μ2=12，標準差σ2=3，相關系數ρ=0.8.確定隨機變量x1，x2的聯合累積分布函數.

1) 選取Clayton Copula，Frank Copula，Gumbel Copula函數構造變量x1,x2的聯合累積分布函數，可以寫成

F(x1,x2)=C(u1,u2;θ).

(15)

2) 通過樣本數據，采取半參數估計法(CML估計)確定參數θ，求得參數θ值如表3所示.

表3 參數θ值

3)由式(3)可以得到x1,x2的聯合累積分布函數為：

在此基礎上，利用文中所提的構造多參數Copula函數方法確定聯合累積分布函數.

同樣利用上述得出的100組數據點，由2.2.1確定擬合樣本，選用OriginLab軟件中非線性曲線擬合，算法選擇正交距離算法，設置算法權重系數均為0.5，參數初始值設為1，迭代一次得到相關參數θ1、θ2如表4所示.

表4 參數θ1、θ2的值

由表4可以得到x1,x2的聯合累積分布函數為：

結果對比

選取顯著性指標α=0.05，對上述6種函數進行K-S擬合檢驗，并利用OLS準則、AIC準則進行函數優選.結果對比如表5.

表5 結果對比

由表可知，文中所構造的3個函數均通過顯著性水平檢驗，擬合效果也優于常用的3類阿基米德Copula函數.其中，③的準則計算值最小，所以③為最優聯合累積分布函數.其表達式如前所示:

由于邊緣分布為正態分布，為驗證已有樣本外數據的準確性，所以采用二維正態分布模型計算所得的理論值與上述6種函數計算結果做對比.隨機生成一組數據，將二維正態分布函數理論值與上述6種函數所計算的值繪制成散點圖，如圖1.

圖1 二維正態分布理論值與6種函數結果對比圖

圖1所示，文中構造的函數①、②、③所繪制的散點均勻地分布在45°對角線附近，說明所構造的函數模型是合理的，也能更直觀的看出其優于傳統的阿基米德Copula函數.

算例2

現有沙潁河流域河南段干流下游周口水文站水文干旱特征變量統計結果，如表6.試根據表中數據構建水文干旱特征變量的二維聯合分布模型.

根據表中數據，選取Frank Copula，Clayton Copula，Gumbel Copula以及所構造的3種函數①、②、③.其中，使用Kolmigrov-Smimov(K-S)來檢驗各聯合累積分布函數模型，運用離差平方和準則(OLS)、赤池信息準則(AIC)對Copula聯合累積分布函數模型進行優選.結果見表7.

表6 周口水文站水文干旱特征變量

續表6

表7 K-S檢驗結果

由表7可知，在顯著性指標α=0.01的情況下，本文構造的聯合累積分布函數①、②、③在D-S、D-P、S-P 3種情況下均通過了K-S檢驗，并且在OLS準則、AIC準則下均優于傳統阿基米德Copula函數.

4 結語

Copula函數理論為構建相關性隨機變量的聯合累積分布函數提供了一種全新的途徑，在已有理論的基礎上，通過乘積生成元構造出了一種新的隱式Copula函數，實現了優于傳統阿基米德Copula方法的擬合效果.針對隱函數中參數求解困難的問題，給出了一種基于經驗分布函數和ODR擬合算法相結合的求解未知參數方法，從而實現了對聯合累積分布函數的構建.最后，通過算例驗證了本文所提方法的有效性，在K-S檢驗以及OLS準則、AIC準則中的表現優于傳統阿基米德Copula函數，實現了對傳統阿基米德Copula函數擬合精度的提高，由于是直接利用樣本的經驗分布函數，因此函數有更大的適用性，并且能隨著樣本數量的提高而自動改良函數的精度.