一種求解相關(guān)隨機(jī)變量聯(lián)合累積分布函數(shù)的方法

楊 永，李海濱

(內(nèi)蒙古工業(yè)大學(xué) 理學(xué)院，內(nèi)蒙古自治區(qū) 呼和浩特 010051)

聯(lián)合概率累積分布函數(shù)是概率論中重要的概念，在工程結(jié)構(gòu)可靠性分析中有著廣泛應(yīng)用.實(shí)際應(yīng)用中，根據(jù)已有的一維函數(shù)模型，如正態(tài)分布、對(duì)數(shù)正態(tài)分布、極值Ι型分布、指數(shù)分布、威布爾分布等，能夠容易得到邊緣累積分布函數(shù)，對(duì)具有獨(dú)立性的數(shù)據(jù)，聯(lián)合累積分布函數(shù)可以通過(guò)邊緣累積分布函數(shù)的直接相乘得到.但在工程中，所測(cè)得的數(shù)據(jù)往往存在相關(guān)性，如巖土工程抗剪強(qiáng)度參數(shù)中的黏聚力和內(nèi)摩擦角、基樁荷載—位移雙曲線等.為了能更好的進(jìn)行可靠性分析，必須考慮數(shù)據(jù)間的相關(guān)性.

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

1959年，Sklar指出，任意一個(gè)多維聯(lián)合累積分布函數(shù)都可以分解為對(duì)應(yīng)的邊緣累積分布函數(shù)和一個(gè)Copula函數(shù)[4].Copula理論的提出為聯(lián)合累積分布函數(shù)的構(gòu)造提供了一個(gè)全新的途徑，由于分開(kāi)構(gòu)造邊緣累積分布函數(shù)和Copula函數(shù)，因此可以構(gòu)造出任意邊緣累積分布函數(shù)和任意相關(guān)結(jié)構(gòu)的聯(lián)合累積分布函數(shù)，而且對(duì)不同相關(guān)結(jié)構(gòu)的數(shù)據(jù)，有多種Copula函數(shù)可以選擇，如Gaussian Copula、t-Copula、Frank Copula等.由于Copula函數(shù)的優(yōu)良性質(zhì)，其在金融、水文、巖土等多個(gè)領(lǐng)域得到了廣泛應(yīng)用[5].而阿基米德Copula函數(shù)是所有Copula函數(shù)中最常用的一類(lèi)函數(shù)，具有形式簡(jiǎn)單、對(duì)稱(chēng)性、可結(jié)合性等優(yōu)點(diǎn)，而且只要找到“生成元”就可以構(gòu)造出這一類(lèi)函數(shù).

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

首先介紹了最常用的阿基米德Copula函數(shù)及其構(gòu)造具有相關(guān)性的二維隨機(jī)變量聯(lián)合累積分布函數(shù)的方法，然后提出一種多參數(shù)Copula函數(shù)的構(gòu)造方法以及對(duì)其構(gòu)造的聯(lián)合累積分布函數(shù)進(jìn)行擬合和參數(shù)估計(jì)，最后利用算例對(duì)其有效性做出驗(yàn)證.

1 基于Copula理論的聯(lián)合概率分布函數(shù)構(gòu)造方法

Copula函數(shù)是對(duì)多元隨機(jī)變量相關(guān)結(jié)構(gòu)的一種刻畫(huà)，可以理解成多元隨機(jī)變量的聯(lián)合累積分布函數(shù)和邊緣累積分布函數(shù)之間聯(lián)系的紐帶[9].

1.1 Copula函數(shù)的定義

Copula理論最早由Sklar于1959年提出.Sklar指出，任意一個(gè)多維聯(lián)合累積分布函數(shù)都可以分解為相應(yīng)的邊緣累積分布函數(shù)和一個(gè)Copula函數(shù)，該Copula函數(shù)確定了變量間的相關(guān)性，包括相關(guān)系數(shù)的大小和相關(guān)結(jié)構(gòu)的類(lèi)型.對(duì)于n維情形,根據(jù)Sklar定理，可將變量x1,x2,…,xn的聯(lián)合累積分布函數(shù)F(x1,x2,…,xn)表示為：

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

(1)

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

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

(2)

(3)

(4)

1.2 Copula函數(shù)的相關(guān)參數(shù)θ

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

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

(5)

其對(duì)數(shù)似然函數(shù)為:

(6)

對(duì)于(6)中的邊緣累積分布函數(shù)值，采用經(jīng)驗(yàn)分布函數(shù)值Femp代替，然后根據(jù)最大似然估計(jì)方法，求解出對(duì)數(shù)似然函數(shù)的最大值點(diǎn)，即可估計(jì)出相關(guān)參數(shù)θ.

(7)

1.3 Copula函數(shù)模型的擬合優(yōu)度檢驗(yàn)

判定選用的函數(shù)模型是否能夠較好的描述變量之間的關(guān)系，需要對(duì)Copula函數(shù)進(jìn)行擬合優(yōu)度檢驗(yàn)，理論上，統(tǒng)計(jì)學(xué)常用的變量分布假設(shè)檢驗(yàn)方法均適用于Copula函數(shù)檢驗(yàn)，如皮爾遜擬合優(yōu)度χ2檢驗(yàn)、Kolmigrov檢驗(yàn)、正態(tài)W檢驗(yàn)等[13].

文中采用Kolmigrov-Smimov(K-S)來(lái)檢驗(yàn)Copula聯(lián)合累積分布函數(shù)模型，用離差平方和準(zhǔn)則(OLS)、赤池信息準(zhǔn)則(AIC)對(duì)Copula聯(lián)合累積分布函數(shù)模型進(jìn)行優(yōu)選.

以二維為例，K-S檢驗(yàn)統(tǒng)計(jì)量D、離差平方和準(zhǔn)則(OLS)、赤池信息準(zhǔn)則(AIC)定義如下：

1) K-S檢驗(yàn)統(tǒng)計(jì)量D

(8)

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

2) 離差平方和準(zhǔn)則(OLS)

(9)

其中，

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

Femp(x1i,x2i)為經(jīng)驗(yàn)分布函數(shù)值，C(u1i,u2i)為理論聯(lián)合累積分布函數(shù)值.

3) 赤池信息準(zhǔn)則(AIC)

(10)

AIC=nln(MSE)+2k,

其中，F(xiàn)emp(x1i,x2i)為經(jīng)驗(yàn)分布函數(shù)值，C(u1i,u2i)為理論聯(lián)合累積分布函數(shù)值，k為模型參數(shù)的個(gè)數(shù).

2 一種乘積生成元構(gòu)造的多參數(shù)Copula

不同的Copula函數(shù)用于描述具有不同相關(guān)特性的相關(guān)變量，在Copula函數(shù)中，阿基米德Copula函數(shù)是現(xiàn)在最常用也是最重要的一類(lèi)Copula，由Genest和Mackey于1986年所提出，其函數(shù)表達(dá)式為：

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

(11)

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

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

表1 阿基米德Copula函數(shù)及其生成元

2.1 一種多參數(shù)Copula函數(shù)的構(gòu)造

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

表2 乘積生成元及其隱式Copula函數(shù)

2.2 相關(guān)參數(shù)θ1、θ2的確定方法

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

2.2.1 擬合數(shù)據(jù)

在統(tǒng)計(jì)學(xué)中，經(jīng)驗(yàn)分布函數(shù)是對(duì)已有樣本累積分布函數(shù)的一種估計(jì)，根據(jù)Glivenko-Cantelli定理，隨著樣本數(shù)的增加，經(jīng)驗(yàn)分布逐步收斂于真實(shí)累積分布值，因此采用邊緣分布的經(jīng)驗(yàn)分布函數(shù)作為邊緣累積分布函數(shù).而對(duì)各個(gè)樣本點(diǎn)處的聯(lián)合累積分布函數(shù)值，同樣可由聯(lián)合經(jīng)驗(yàn)分布函數(shù)代替.以前者為自變量，后者為因變量，即可構(gòu)成擬合樣本.以n組二維數(shù)據(jù)樣本(x1,x2)為例：

邊緣累積分布函數(shù)

(12)

其中，

{x11,x12,…，x1n}*表示集合{x11,x12,…，x1n}中不大于x1i的個(gè)數(shù),

{x21,x22,…，x2n}*表示集合{x21,x22,…，x2n}中不大于x2i的個(gè)數(shù).

聯(lián)合累積分布函數(shù)：

(13)

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

2.2.2 正交距離曲線擬合

由于采用了經(jīng)驗(yàn)分布函數(shù)，為考慮邊緣經(jīng)驗(yàn)分布函數(shù)代邊緣累積分布函數(shù)所帶來(lái)的誤差，即自變量誤差，所以文中算法選取正交最小二乘法.正交最小二乘算法又稱(chēng)正交距離曲線擬合，其與普通最小二乘的區(qū)別是考慮了函數(shù)擬合過(guò)程中自變量的誤差，以正交距離的殘差平方和極小為準(zhǔn)則進(jìn)行曲線擬合，使擬合結(jié)果從整體上達(dá)到擬合最佳[14].

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

(14)

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

3 數(shù)值算例

算例1

隨機(jī)生成100組數(shù)據(jù)點(diǎn)作為樣本.其中，隨機(jī)變量x1為正態(tài)分布，均值μ1=8，標(biāo)準(zhǔn)差σ1=2，隨機(jī)變量x2也為正態(tài)分布，均值μ2=12，標(biāo)準(zhǔn)差σ2=3，相關(guān)系數(shù)ρ=0.8.確定隨機(jī)變量x1，x2的聯(lián)合累積分布函數(shù).

1) 選取Clayton Copula，F(xiàn)rank Copula，Gumbel Copula函數(shù)構(gòu)造變量x1,x2的聯(lián)合累積分布函數(shù)，可以寫(xiě)成

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

(15)

2) 通過(guò)樣本數(shù)據(jù)，采取半?yún)?shù)估計(jì)法(CML估計(jì))確定參數(shù)θ，求得參數(shù)θ值如表3所示.

表3 參數(shù)θ值

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

在此基礎(chǔ)上，利用文中所提的構(gòu)造多參數(shù)Copula函數(shù)方法確定聯(lián)合累積分布函數(shù).

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

表4 參數(shù)θ1、θ2的值

由表4可以得到x1,x2的聯(lián)合累積分布函數(shù)為：

結(jié)果對(duì)比

選取顯著性指標(biāo)α=0.05，對(duì)上述6種函數(shù)進(jìn)行K-S擬合檢驗(yàn)，并利用OLS準(zhǔn)則、AIC準(zhǔn)則進(jìn)行函數(shù)優(yōu)選.結(jié)果對(duì)比如表5.

表5 結(jié)果對(duì)比

由表可知，文中所構(gòu)造的3個(gè)函數(shù)均通過(guò)顯著性水平檢驗(yàn)，擬合效果也優(yōu)于常用的3類(lèi)阿基米德Copula函數(shù).其中，③的準(zhǔn)則計(jì)算值最小，所以③為最優(yōu)聯(lián)合累積分布函數(shù).其表達(dá)式如前所示:

由于邊緣分布為正態(tài)分布，為驗(yàn)證已有樣本外數(shù)據(jù)的準(zhǔn)確性，所以采用二維正態(tài)分布模型計(jì)算所得的理論值與上述6種函數(shù)計(jì)算結(jié)果做對(duì)比.隨機(jī)生成一組數(shù)據(jù)，將二維正態(tài)分布函數(shù)理論值與上述6種函數(shù)所計(jì)算的值繪制成散點(diǎn)圖，如圖1.

圖1 二維正態(tài)分布理論值與6種函數(shù)結(jié)果對(duì)比圖

圖1所示，文中構(gòu)造的函數(shù)①、②、③所繪制的散點(diǎn)均勻地分布在45°對(duì)角線附近，說(shuō)明所構(gòu)造的函數(shù)模型是合理的，也能更直觀的看出其優(yōu)于傳統(tǒng)的阿基米德Copula函數(shù).

算例2

現(xiàn)有沙潁河流域河南段干流下游周口水文站水文干旱特征變量統(tǒng)計(jì)結(jié)果，如表6.試根據(jù)表中數(shù)據(jù)構(gòu)建水文干旱特征變量的二維聯(lián)合分布模型.

根據(jù)表中數(shù)據(jù)，選取Frank Copula，Clayton Copula，Gumbel Copula以及所構(gòu)造的3種函數(shù)①、②、③.其中，使用Kolmigrov-Smimov(K-S)來(lái)檢驗(yàn)各聯(lián)合累積分布函數(shù)模型，運(yùn)用離差平方和準(zhǔn)則(OLS)、赤池信息準(zhǔn)則(AIC)對(duì)Copula聯(lián)合累積分布函數(shù)模型進(jìn)行優(yōu)選.結(jié)果見(jiàn)表7.

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

續(xù)表6

表7 K-S檢驗(yàn)結(jié)果

由表7可知，在顯著性指標(biāo)α=0.01的情況下，本文構(gòu)造的聯(lián)合累積分布函數(shù)①、②、③在D-S、D-P、S-P 3種情況下均通過(guò)了K-S檢驗(yàn)，并且在OLS準(zhǔn)則、AIC準(zhǔn)則下均優(yōu)于傳統(tǒng)阿基米德Copula函數(shù).

4 結(jié)語(yǔ)

Copula函數(shù)理論為構(gòu)建相關(guān)性隨機(jī)變量的聯(lián)合累積分布函數(shù)提供了一種全新的途徑，在已有理論的基礎(chǔ)上，通過(guò)乘積生成元構(gòu)造出了一種新的隱式Copula函數(shù)，實(shí)現(xiàn)了優(yōu)于傳統(tǒng)阿基米德Copula方法的擬合效果.針對(duì)隱函數(shù)中參數(shù)求解困難的問(wèn)題，給出了一種基于經(jīng)驗(yàn)分布函數(shù)和ODR擬合算法相結(jié)合的求解未知參數(shù)方法，從而實(shí)現(xiàn)了對(duì)聯(lián)合累積分布函數(shù)的構(gòu)建.最后，通過(guò)算例驗(yàn)證了本文所提方法的有效性，在K-S檢驗(yàn)以及OLS準(zhǔn)則、AIC準(zhǔn)則中的表現(xiàn)優(yōu)于傳統(tǒng)阿基米德Copula函數(shù)，實(shí)現(xiàn)了對(duì)傳統(tǒng)阿基米德Copula函數(shù)擬合精度的提高，由于是直接利用樣本的經(jīng)驗(yàn)分布函數(shù)，因此函數(shù)有更大的適用性，并且能隨著樣本數(shù)量的提高而自動(dòng)改良函數(shù)的精度.