牛巖溪，梁馮珍

（天津大學(xué) 數(shù)學(xué)學(xué)院，天津 300350）

0 引言

貝葉斯網(wǎng)絡(luò)(Bayesian networks)也稱信念網(wǎng)絡(luò)(Belief networks)，是描述隨機變量間依賴關(guān)系的圖形模式，被廣泛用于不確定性問題的智能化求解。該網(wǎng)絡(luò)由有向無環(huán)圖(DAG)和概率參數(shù)構(gòu)成。DAG的節(jié)點表示隨機變量，節(jié)點間的有向弧表示變量間的相依關(guān)系，具體的相依程度則通過概率參數(shù)反映。1988年，Pearl[1]首先給出其嚴格定義，并創(chuàng)建貝葉斯網(wǎng)絡(luò)基礎(chǔ)理論體系。貝葉斯網(wǎng)絡(luò)具有多功能性、有效性和開放性等特征，能夠有效地轉(zhuǎn)化數(shù)據(jù)為知識，并利用這些知識進行推理，以解決分析、預(yù)測和控制等方面的問題。Kurowicka和Cooke(2002)[2]將Pair-Copula和貝葉斯網(wǎng)絡(luò)結(jié)合，在繼承了一般貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)的優(yōu)點的同時，結(jié)合Copula函數(shù)的優(yōu)良性質(zhì)，建立一種新的多元統(tǒng)計模型——Pair-Copula貝葉斯網(wǎng)絡(luò)(PCBNs)模型。PCBNs模型通過結(jié)構(gòu)推斷和相關(guān)性分析，反映變量之間的因果相依關(guān)系，為各領(lǐng)域研究人員和管理者提供決策支持，近年來已在金融、工程等領(lǐng)域有所應(yīng)用。

構(gòu)建Pair-Copula貝葉斯網(wǎng)絡(luò)模型包括結(jié)構(gòu)估計和參數(shù)估計。Bauer等(2012)[3]給出了DAG模型下Pair-Copula的ML參數(shù)估計的方法。結(jié)構(gòu)估計常用的機器學(xué)習(xí)方法有PC算法[4]、爬山(Hill-climbing)算法等。Colombo(2012)[5]對經(jīng)典PC算法模型進行了改進，降低了檢驗順序?qū)z驗結(jié)果的影響。然而，在結(jié)構(gòu)估計中，核心步驟是變量之間條件獨立性的檢驗。條件獨立性檢驗的難度要遠遠大于非條件獨立性檢驗。Bouezmarni等（2009）[6]基于Hellinger距離，結(jié)合Bernstein Copula構(gòu)造了檢驗統(tǒng)計量來檢驗條件獨立性。Zhang（2011）[7]提出基于核函數(shù)的條件獨立檢驗方法。Ramsey等（2014）[8]對核函數(shù)模型進行了簡化和推廣。

本文在改進的PC算法的基礎(chǔ)上，提出基于多項式回歸殘差的條件獨立性檢驗方法，并進行仿真模擬實驗以及氣象數(shù)據(jù)實證分析。該方法可以良好地檢驗變量間的條件獨立關(guān)系，使用網(wǎng)絡(luò)結(jié)構(gòu)體現(xiàn)變量間的相依和獨立關(guān)系，并結(jié)合Pair-Copula得到完整的相依關(guān)系推斷模型以及相應(yīng)的聯(lián)合密度函數(shù)。

1 Pair-Copula貝葉斯網(wǎng)絡(luò)模型

當變量集給定時，構(gòu)建Pair-Copula貝葉斯網(wǎng)絡(luò)主要包括結(jié)構(gòu)估計和參數(shù)估計。Pair-Copula貝葉斯網(wǎng)絡(luò)允許在變量的子集間定義條件獨立性，這也是貝葉斯網(wǎng)絡(luò)最重要的性質(zhì)之一。在結(jié)構(gòu)估計中，該性質(zhì)將變量間的條件獨立關(guān)系從復(fù)雜的網(wǎng)絡(luò)中除去，大大簡化了網(wǎng)絡(luò)模型結(jié)構(gòu)和相應(yīng)的聯(lián)合密度函數(shù)的復(fù)雜程度，使得變量間的相依和獨立關(guān)系可以通過簡潔的有向無環(huán)圖直觀地展現(xiàn)。檢驗變量間的條件獨立性在Pair-Copula貝葉斯網(wǎng)絡(luò)的構(gòu)建和相依關(guān)系的發(fā)現(xiàn)中尤為重要。

1.1 馬爾科夫概率測度

在圖形模式中，和變量條件獨立性相對應(yīng)的是節(jié)點的D-馬爾科夫?qū)傩?。令V≠?為一個有限集合，E?E={(v,w)∈V×V|v≠w}。則D=(V,E)表示一個由頂點集V和邊集E構(gòu)成的DAG，Dm為D的道德圖[3]。對于任意v∈V，令：

令P為?d上的概率測度，d=|V|。X為?d上概率分布為P的d維隨機變量。對于I?V，記XI=(Xv)v∈I。若對于兩兩互不相交的集合I,J,K?V，有XI和XJ在給定XK的條件下獨立，則記為XI⊥XJ|XK。

當對于?v∈V，有：

則稱P具有局部D-馬爾科夫?qū)傩浴?/p>

當對于兩兩互不相交的集合I,J,K?V，有：

則稱P具有全局D-馬爾科夫?qū)傩?，其中An(?)表示最小祖先集。當概率測度滿足式（1）和式（2）時，則稱P具有D-馬爾科夫?qū)傩浴-馬爾科夫?qū)傩钥梢愿M一步地展示出由有向無環(huán)圖D表示的條件獨立性。

現(xiàn)令P有Lebesgue概率密度f，則P是D-馬爾科夫的當且僅當f有如下形式的D-遞歸因式分解：

其中fv|pa(v)(?|xpa(v))表示在給定Xpa(v)=xpa(v)的條件下，Xv的條件概率密度函數(shù)[3]。

1.2 Pair-Copula貝葉斯網(wǎng)絡(luò)模型

令D=(V,E)為一個DAG，P為?d上絕對連續(xù)的D-馬爾科夫概率測度，d=|V|。X為?d上概率分布為P的d維隨機變量。對于?v∈V，s.t.|pa(v)|≥1 ，令wv:{1,…,|，是一個雙射。令＜v為pa(v)上的一個全序關(guān)系，在這個關(guān)系中，?i,j∈{1,…,有i＜j當且僅當wi＜vwj。則?v∈V，w∈pa(v)，記：

由Sklar定理[9]可知，P的概率密度函數(shù)可以唯一地分解成一元邊緣密度函數(shù)fi和Copula函數(shù)c的形式。Bauer等(2012)[3]指出c又可以進一步分解成（條件）成對Copulacv,w|pa(v;w)的形式，其中每個條件Copula關(guān)系對應(yīng)于DAG中的一條邊w→v，cv,w|pa(v;w)表示相應(yīng)密度函數(shù)。因此，P的概率密度函數(shù)f可以寫成如下形式：

其中x=(xv)v∈V∈?d。

令u=(u1,…,un)，n?? ，是 [0,1]d上的隨機變量U的i.i.d.觀測值，其Copula分布族為{CD;θ|θ∈Θ}，邊緣分布為均勻分布。則式（3）對應(yīng)的對數(shù)似然函數(shù)為：

同樣作為由Pair-Copula構(gòu)成的圖結(jié)構(gòu)，Vine-Copula結(jié)構(gòu)也是數(shù)據(jù)建模中常用模型之一。Vine-Copula方法隨著變量維度的增加，可選結(jié)構(gòu)的種類以及待估參數(shù)數(shù)量將隨之以平方函數(shù)速度增加，運算量較大。從式（3）的聯(lián)合密度函數(shù)表達式可以看出，Pair-Copula貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)與Vine-Copula結(jié)構(gòu)相似，但卻剔除了條件獨立變量之間的關(guān)系。易知，建立一個d維Regular vine(R-vine)模型結(jié)構(gòu)[9]，需要定義的Pair-Copula數(shù)量為而由于考慮了條件獨立性，Pair-Copula貝葉斯網(wǎng)絡(luò)結(jié)構(gòu)所需Pair-Copula數(shù)量將減少。隨著維度增加，這個優(yōu)勢會更加凸顯。

1.3 基于多項式回歸的條件獨立性檢驗

給定變量集，構(gòu)建Pair-Copula貝葉斯網(wǎng)絡(luò)的首要任務(wù)即利用變量間的馬爾科夫?qū)傩宰R別DAGD=(V,E)。因此，檢驗變量間的條件獨立性在貝葉斯網(wǎng)絡(luò)的構(gòu)建和相依關(guān)系的發(fā)現(xiàn)中尤為重要。改進的PC算法[5]是一種數(shù)據(jù)驅(qū)動下的網(wǎng)絡(luò)結(jié)構(gòu)估計算法。在此算法的基礎(chǔ)上，本文提出一種基于多項式回歸判斷變量間條件獨立性的方法。

根據(jù)數(shù)學(xué)分析中的知識易知，任何函數(shù)都可以近似地用多項式表示。因此，兩個變量間的關(guān)系可以用多項式進行逼近，即多項式回歸，變量無關(guān)部分則通過回歸殘差反映。相較于核回歸，多項式回歸在處理實際問題中更為常用，使用方便且容易解釋。多項式回歸模型的一般形式為：

下面給出在多項式回歸模型下的條件獨立性判斷方法。

設(shè)X,Y,Z為?上的隨機變量，X和Y在給定Z的條件下獨立，當且僅當?shù)呔S時，估計密度函數(shù)困難。因此，本文避開明確的密度估計，考慮X,Y,Z間的非線性回歸，而不做其他的分布或函數(shù)結(jié)構(gòu)的假設(shè)。

本文選用多項式回歸模型，認為X和Y在給定Z的條件下獨立，當且僅當Z回歸到X的殘差與Z回歸到Y(jié)的殘差之間相互獨立。即：

其中SX和SY為多項式函數(shù)，則：

不同于非條件獨立關(guān)系，當給定Z時，X和Y的關(guān)系無法直接獲得。通過多項式回歸，分別從X和Y中將Z的影響分離，與Z不相關(guān)的部分放在誤差變量ε中。此時，本文將復(fù)雜條件獨立性的判斷轉(zhuǎn)化為了非條件獨立性的判斷。不同于數(shù)理統(tǒng)計中關(guān)于兩變量獨立的嚴格定義，在實際應(yīng)用中，絕對獨立的情況并不常見，因此在本文所建模型之下，認為給定Z時，X和Y不相關(guān)即可判定為條件獨立。此處可以選取任意一種合理的相關(guān)系數(shù)來檢驗ε1,ε2之間是否相關(guān)。

1.4 基于多項式回歸的PCBNs模型構(gòu)建

獲取實際數(shù)據(jù)，對數(shù)據(jù)做一定的預(yù)處理，構(gòu)建Pair-Copula貝葉斯網(wǎng)絡(luò)，主要步驟如下：

步驟1：結(jié)構(gòu)估計。

對樣本使用改進PC算法進行歷遍，其中獨立性檢驗步驟改為上述多項式回歸方法。由所有節(jié)點的完全無向圖[3]出發(fā)，假設(shè)節(jié)點順序未知。

若εiK⊥εjK，接受原假設(shè)H0，刪去邊i-j；否則拒絕原假設(shè)H0，保留邊i-j。此時，條件獨立性的判斷被轉(zhuǎn)化為相對簡單的非條件獨立性的判斷，本文選取Kendall相關(guān)系數(shù)(Kendall’sτ)，若|τ|較小，εiK與εjK不相關(guān)，則判斷

此時DAG對應(yīng)的無向圖[3]得到，記為DU。

步驟2：方向確定。

無向圖DU確定后，需確定邊的方向。設(shè)Sij為i,j的分割集，i≠j∈V，(i,j)?EDU，(i,j)?EDU。

對于i∈V，j?ad(i)，k∈ad(i)∩ad(j)，如果k?Sij，則i-k-j的方向為i→k←j。（v-結(jié)構(gòu)）

當尋找到所有存在的v-結(jié)構(gòu)后，其他無向邊方向規(guī)則為：

若DU中包含i→j且k?ad(i)，則j-k方向為j→k（否則將存在新的v-結(jié)構(gòu)）。

若DU中包含i→k→j，則i-j的方向為i→j（否則將存在環(huán)結(jié)構(gòu)）。

若DU中包含i→k→j和i→I→j，l?ad(k)，則i-j的方向為i→j（否則將存在新的v-結(jié)構(gòu)或環(huán)結(jié)構(gòu)）。

只使用馬爾科夫性質(zhì)，有時不足以得到完整的相依關(guān)系，可能存在某些邊的方向無法確定的情況，但可作為進一步探索的起點。實際應(yīng)用中，可以根據(jù)數(shù)據(jù)的具體情況，結(jié)合其他方法和經(jīng)驗確定或修改方向。所有邊方向確定后，得到有向無環(huán)圖D對應(yīng)的鏈圖[3]D*。

步驟3：參數(shù)估計。

確定邊對應(yīng)的條件相依關(guān)系。結(jié)合Vine-Copula結(jié)構(gòu)的參數(shù)估計方法，為各邊構(gòu)造藤結(jié)構(gòu)，進行參數(shù)估計[3]。本文選取D-vine結(jié)構(gòu)，首先求樣本的經(jīng)驗分布函數(shù)，用于Copula建模。由于網(wǎng)絡(luò)結(jié)構(gòu)的方向確定，因此可以確定每個節(jié)點的祖先集，對于D中i∈V，j∈ad(i)，k?an(i)∩an(j)，根據(jù)D-vine的參數(shù)估計方法，為邊構(gòu)造一個或多個藤結(jié)構(gòu)，估計所有(xi,xj)|xk的Pair-Copula類型及參數(shù)，求得參數(shù)對應(yīng)的Kendall’sτ值，結(jié)合實際問題選取|τ|較大的作為該邊的相關(guān)關(guān)系。由于本文主要使用獨立性關(guān)系作為抽樣依據(jù)，對于相依關(guān)系的依賴較少，因此條件集的選擇范圍也較寬，在實際問題應(yīng)用中十分靈活。Pair-Copula函數(shù)類型和參數(shù)確定后，即可寫出相應(yīng)的聯(lián)合密度函數(shù)，并得到最終的Pair-Copula貝葉斯網(wǎng)絡(luò)。

2 仿真模擬

本文考慮一個簡單的5節(jié)點貝葉斯網(wǎng)絡(luò)，樣本容量為5000，使用R語言編寫程序。假設(shè)給定網(wǎng)絡(luò)模型以及變量間的相依關(guān)系如圖1所示。在該網(wǎng)絡(luò)中，變量間的條件獨立關(guān)系為2⊥3|1，1⊥4|(2,3)，1⊥5|4，2⊥5|4，3⊥5|4 。首先確定抽樣方法，對樣本使用上次中的步驟，構(gòu)建完整Pair-Copula貝葉斯網(wǎng)絡(luò)。

圖1 5節(jié)點貝葉斯網(wǎng)絡(luò)

2.1 抽樣方法

根據(jù)條件獨立關(guān)系和相關(guān)關(guān)系構(gòu)造變量之間的多項式組，用于抽樣。由于任何函數(shù)都可用（分段）多項式形式表示，因此可以任意構(gòu)造能體現(xiàn)相應(yīng)的獨立、相關(guān)關(guān)系多項式。在多項式回歸中，當自變量的冪超過3時，回歸系數(shù)的解釋將變得困難，回歸函數(shù)也變得很不穩(wěn)定，對回歸模型的應(yīng)用會受到影響。因而，冪次超過3的模型不常使用。本文以二次方多項式為例，建立如下多項式關(guān)系：

抽樣步驟如下：

步驟1：根據(jù)殘差的獨立關(guān)系，利用獨立高斯Copula生成殘差ε1,ε2,ε3,ε4。

設(shè)ε1,ε2,ε3,ε4為 回 歸 殘 差 ，ε1⊥ε2，則 (ε1,ε2)～GassianCopula(0,dim=2)，生成兩組聯(lián)合分布為獨立高斯Copula的隨機數(shù)，分別賦予ε1,ε2。生成（0,1）上均勻分布隨機數(shù)w3～U(0,1)，w4～U(0,1)，根據(jù)ε2⊥ε3，有為獨立高斯 Copula的條件逆Copula函數(shù)。同理

步驟2：根據(jù)多項式關(guān)系生成樣本。首先需生成（0,1）上均勻分布隨機數(shù)u1。

此時模擬樣本得到。

2.2 結(jié)構(gòu)估計

根據(jù)1.4中的步驟1，由所有節(jié)點的完全無向圖出發(fā)檢驗條件獨立性，如圖2所示。使用基于多項式回歸的改進PC算法進行節(jié)點歷遍后，得到邊的具體刪留情況如表1所示，此時DAG對應(yīng)的無向圖確定，記為DU。

圖2完全無向圖

表1 完全無向圖邊的刪留情況

在該例中，5節(jié)點結(jié)構(gòu)到此歷遍結(jié)束。對于更高維的數(shù)據(jù)，可以不考慮條件集長度大于3的情況，或僅作為參考，綜合前幾層歷遍決定邊的刪留。

根據(jù)1.4中的步驟2，當通過v-結(jié)構(gòu)及其他無向邊方向規(guī)則為無向圖確定方向后，得到鏈圖D*，如圖3所示，此時仍有兩條邊不能完全確定方向，即1—2和1—3。實際問題中可根據(jù)情況或結(jié)合其他方法確定。暫且假定為1→2，1→3。此時網(wǎng)絡(luò)完整結(jié)構(gòu)得到，如圖4所示。

圖3 鏈圖D*

圖4有向無環(huán)圖D

2.3 參數(shù)估計

網(wǎng)絡(luò)結(jié)構(gòu)確定后，并不能最終確定每條邊對應(yīng)的條件相依關(guān)系。如邊4→5，對應(yīng)的關(guān)系可能是45|12、45|123或其他。根據(jù)1.4中的步驟3，為邊構(gòu)造藤結(jié)構(gòu)確定Copula類型，估計參數(shù)。本文只考慮Copula函數(shù)中最常見的幾種函數(shù)類型。

最終得到的Pair-Copula類型及參數(shù)如表2所示：

表2 Copula類型及參數(shù)

聯(lián)合密度函數(shù)為：

Pair-Copula貝葉斯網(wǎng)絡(luò)得到，如圖5所示：

圖5 Pair-Copula貝葉斯網(wǎng)絡(luò)

3 實證

在空調(diào)系統(tǒng)設(shè)計過程，室外氣象參數(shù)是負荷計算、設(shè)備選型等必需的基礎(chǔ)參數(shù)，直接影響系統(tǒng)設(shè)計和運行效果。氣象參數(shù)之間具有一定程度的相依關(guān)系，本文選取四種典型參數(shù)：干球溫度、濕球溫度、含濕量、太陽輻射，將上文介紹的方法應(yīng)用到逐時氣象數(shù)據(jù)中進行實證分析，探討這四種氣象參數(shù)之間的相關(guān)性。數(shù)據(jù)來自天津市氣候中心，選取天津市2010年全年逐時干球溫度、濕球溫度、含濕量、太陽輻射觀測值，每項觀測值包含8760個數(shù)值。

設(shè)干球溫度、濕球溫度、含濕量、太陽輻射依次為變量x1,x2,x3,x4。根據(jù)物理意義及現(xiàn)實情況，溫度、濕度等空氣狀態(tài)受太陽輻射強度影響，反之，太陽輻射不受溫度等影響，因此確定太陽輻射為首層節(jié)點，不考慮在其他變量的條件下太陽輻射與某一變量的關(guān)系，即只考慮太陽輻射與其他變量的無條件關(guān)系。其他變量的節(jié)點順序未知。

首先使用基于多項式回歸的改進PC算法進行節(jié)點歷遍，選取|τ|=0.3 ，若|τ|＜0.3 認為變量間（條件）不相關(guān)，判定滿足條件獨立性。得到邊的具體刪留情況如表3所示。

表3 氣象數(shù)據(jù)完全無向圖邊的刪留情況

按照方向規(guī)則為無向邊確定方向，得到DAG如下頁圖6所示。

圖6氣象數(shù)據(jù)有向無環(huán)圖

網(wǎng)絡(luò)結(jié)構(gòu)確定后，將數(shù)據(jù)變換為Copula水平。由于本文重在研究相關(guān)結(jié)構(gòu)，不對邊緣分布做具體研究，因此通過經(jīng)驗概率積分變換將逐時干球溫度、濕球溫度、含濕量、太陽輻射變換為Copula數(shù)據(jù)?？紤]每條邊可能對應(yīng)的所有（條件）關(guān)系，估計Pair-Copula類型及參數(shù)，并求得參數(shù)對應(yīng)的Kendall’sτ值，結(jié)合實際情況和τ值，選取該邊的相關(guān)關(guān)系。最終得到的Pair-Copula類型及參數(shù)如表4所示：

表4 氣象數(shù)據(jù)Copula類型及參數(shù)

Pair-Copula貝葉斯網(wǎng)絡(luò)如圖7所示。

圖7氣象數(shù)據(jù)Pair-Copula貝葉斯網(wǎng)絡(luò)

從上述Pair-Copula貝葉斯網(wǎng)絡(luò)看出，干球溫度與含濕量、太陽輻射有關(guān)；太陽輻射給定時，濕球溫度與含濕量相關(guān)；含濕量給定時，干濕球溫度相關(guān)。從物理意義角度分析，干球溫度是溫度計自由地被暴露在空氣中所測量的溫度，通常被視作所測量空氣的實際溫度，與含濕量、太陽輻射有關(guān)。濕球溫度是溫度計的球體表面附著有水時，水份蒸發(fā)帶走熱量后球體的溫度。太陽輻射給定時，水的蒸發(fā)量跟空氣的濕度有關(guān)，空氣濕度越大蒸發(fā)量越小，帶走的熱量越少，干濕球溫度差異越??；空氣濕度越小水蒸發(fā)量越大，帶走的熱量也越大，干濕球溫差也就越大。通過該Pair-Copula貝葉斯網(wǎng)絡(luò)可以推斷當某一變量發(fā)生變化時，干球溫度的變化情況。

4 總結(jié)

Pair-Copula貝葉斯網(wǎng)絡(luò)能夠有效地轉(zhuǎn)化數(shù)據(jù)為知識，并利用這些知識進行推理，以解決因果推斷等問題。在改進的PC算法的基礎(chǔ)上，多項式回歸殘差的條件獨立性檢驗方法可以良好地檢驗變量間的條件獨立關(guān)系，結(jié)合Pair-Copula得到完整的相依關(guān)系推斷模型及聯(lián)合密度函數(shù)。在節(jié)點個數(shù)較少的情況下，該方法簡單有效，適用于在實際應(yīng)用中解決金融、工程等領(lǐng)域的相關(guān)問題。當節(jié)點個數(shù)較多時，可綜合其他方法進一步研究探討。