基于多项式回归的Pair-Copula贝叶斯网络模型

牛岩溪，梁冯珍

（天津大学 数学学院，天津 300350）

0 引言

贝叶斯网络(Bayesian networks)也称信念网络(Belief networks)，是描述随机变量间依赖关系的图形模式，被广泛用于不确定性问题的智能化求解。该网络由有向无环图(DAG)和概率参数构成。DAG的节点表示随机变量，节点间的有向弧表示变量间的相依关系，具体的相依程度则通过概率参数反映。1988年，Pearl[1]首先给出其严格定义，并创建贝叶斯网络基础理论体系。贝叶斯网络具有多功能性、有效性和开放性等特征，能够有效地转化数据为知识，并利用这些知识进行推理，以解决分析、预测和控制等方面的问题。Kurowicka和Cooke(2002)[2]将Pair-Copula和贝叶斯网络结合，在继承了一般贝叶斯网络结构的优点的同时，结合Copula函数的优良性质，建立一种新的多元统计模型——Pair-Copula贝叶斯网络(PCBNs)模型。PCBNs模型通过结构推断和相关性分析，反映变量之间的因果相依关系，为各领域研究人员和管理者提供决策支持，近年来已在金融、工程等领域有所应用。

构建Pair-Copula贝叶斯网络模型包括结构估计和参数估计。Bauer等(2012)[3]给出了DAG模型下Pair-Copula的ML参数估计的方法。结构估计常用的机器学习方法有PC算法[4]、爬山(Hill-climbing)算法等。Colombo(2012)[5]对经典PC算法模型进行了改进，降低了检验顺序对检验结果的影响。然而，在结构估计中，核心步骤是变量之间条件独立性的检验。条件独立性检验的难度要远远大于非条件独立性检验。Bouezmarni等（2009）[6]基于Hellinger距离，结合Bernstein Copula构造了检验统计量来检验条件独立性。Zhang（2011）[7]提出基于核函数的条件独立检验方法。Ramsey等（2014）[8]对核函数模型进行了简化和推广。

本文在改进的PC算法的基础上，提出基于多项式回归残差的条件独立性检验方法，并进行仿真模拟实验以及气象数据实证分析。该方法可以良好地检验变量间的条件独立关系，使用网络结构体现变量间的相依和独立关系，并结合Pair-Copula得到完整的相依关系推断模型以及相应的联合密度函数。

1 Pair-Copula贝叶斯网络模型

当变量集给定时，构建Pair-Copula贝叶斯网络主要包括结构估计和参数估计。Pair-Copula贝叶斯网络允许在变量的子集间定义条件独立性，这也是贝叶斯网络最重要的性质之一。在结构估计中，该性质将变量间的条件独立关系从复杂的网络中除去，大大简化了网络模型结构和相应的联合密度函数的复杂程度，使得变量间的相依和独立关系可以通过简洁的有向无环图直观地展现。检验变量间的条件独立性在Pair-Copula贝叶斯网络的构建和相依关系的发现中尤为重要。

1.1 马尔科夫概率测度

在图形模式中，和变量条件独立性相对应的是节点的D-马尔科夫属性。令V≠∅为一个有限集合，E⊆E={(v,w)∈V×V|v≠w}。则D=(V,E)表示一个由顶点集V和边集E构成的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-马尔科夫属性。

当对于两两互不相交的集合I,J,K⊆V，有：

则称P具有全局D-马尔科夫属性，其中An(⋅)表示最小祖先集。当概率测度满足式（1）和式（2）时，则称P具有D-马尔科夫属性。D-马尔科夫属性可以更进一步地展示出由有向无环图D表示的条件独立性。

现令P有Lebesgue概率密度f，则P是D-马尔科夫的当且仅当f有如下形式的D-递归因式分解：

其中fv|pa(v)(⋅|xpa(v))表示在给定Xpa(v)=xpa(v)的条件下，Xv的条件概率密度函数[3]。

1.2 Pair-Copula贝叶斯网络模型

令D=(V,E)为一个DAG，P为ℝd上绝对连续的D-马尔科夫概率测度，d=|V|。X为ℝd上概率分布为P的d维随机变量。对于∀v∈V，s.t.|pa(v)|≥1 ，令wv:{1,…,|，是一个双射。令＜v为pa(v)上的一个全序关系，在这个关系中，∀i,j∈{1,…,有i＜j当且仅当wi＜vwj。则∀v∈V，w∈pa(v)，记：

由Sklar定理[9]可知，P的概率密度函数可以唯一地分解成一元边缘密度函数fi和Copula函数c的形式。Bauer等(2012)[3]指出c又可以进一步分解成（条件）成对Copulacv,w|pa(v;w)的形式，其中每个条件Copula关系对应于DAG中的一条边w→v，cv,w|pa(v;w)表示相应密度函数。因此，P的概率密度函数f可以写成如下形式：

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

令u=(u1,…,un)，nϵℕ ，是 [0,1]d上的随机变量U的i.i.d.观测值，其Copula分布族为{CD;θ|θ∈Θ}，边缘分布为均匀分布。则式（3）对应的对数似然函数为：

同样作为由Pair-Copula构成的图结构，Vine-Copula结构也是数据建模中常用模型之一。Vine-Copula方法随着变量维度的增加，可选结构的种类以及待估参数数量将随之以平方函数速度增加，运算量较大。从式（3）的联合密度函数表达式可以看出，Pair-Copula贝叶斯网络结构与Vine-Copula结构相似，但却剔除了条件独立变量之间的关系。易知，建立一个d维Regular vine(R-vine)模型结构[9]，需要定义的Pair-Copula数量为而由于考虑了条件独立性，Pair-Copula贝叶斯网络结构所需Pair-Copula数量将减少。随着维度增加，这个优势会更加凸显。

1.3 基于多项式回归的条件独立性检验

给定变量集，构建Pair-Copula贝叶斯网络的首要任务即利用变量间的马尔科夫属性识别DAGD=(V,E)。因此，检验变量间的条件独立性在贝叶斯网络的构建和相依关系的发现中尤为重要。改进的PC算法[5]是一种数据驱动下的网络结构估计算法。在此算法的基础上，本文提出一种基于多项式回归判断变量间条件独立性的方法。

根据数学分析中的知识易知，任何函数都可以近似地用多项式表示。因此，两个变量间的关系可以用多项式进行逼近，即多项式回归，变量无关部分则通过回归残差反映。相较于核回归，多项式回归在处理实际问题中更为常用，使用方便且容易解释。多项式回归模型的一般形式为：

下面给出在多项式回归模型下的条件独立性判断方法。

设X,Y,Z为ℝ上的随机变量，X和Y在给定Z的条件下独立，当且仅当但高维时，估计密度函数困难。因此，本文避开明确的密度估计，考虑X,Y,Z间的非线性回归，而不做其他的分布或函数结构的假设。

本文选用多项式回归模型，认为X和Y在给定Z的条件下独立，当且仅当Z回归到X的残差与Z回归到Y的残差之间相互独立。即：

其中SX和SY为多项式函数，则：

不同于非条件独立关系，当给定Z时，X和Y的关系无法直接获得。通过多项式回归，分别从X和Y中将Z的影响分离，与Z不相关的部分放在误差变量ε中。此时，本文将复杂条件独立性的判断转化为了非条件独立性的判断。不同于数理统计中关于两变量独立的严格定义，在实际应用中，绝对独立的情况并不常见，因此在本文所建模型之下，认为给定Z时，X和Y不相关即可判定为条件独立。此处可以选取任意一种合理的相关系数来检验ε1,ε2之间是否相关。

1.4 基于多项式回归的PCBNs模型构建

获取实际数据，对数据做一定的预处理，构建Pair-Copula贝叶斯网络，主要步骤如下：

步骤1：结构估计。

对样本使用改进PC算法进行历遍，其中独立性检验步骤改为上述多项式回归方法。由所有节点的完全无向图[3]出发，假设节点顺序未知。

若εiK⊥εjK，接受原假设H0，删去边i-j；否则拒绝原假设H0，保留边i-j。此时，条件独立性的判断被转化为相对简单的非条件独立性的判断，本文选取Kendall相关系数(Kendall’sτ)，若|τ|较小，εiK与εjK不相关，则判断

此时DAG对应的无向图[3]得到，记为DU。

步骤2：方向确定。

无向图DU确定后，需确定边的方向。设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-结构）

当寻找到所有存在的v-结构后，其他无向边方向规则为：

若DU中包含i→j且k∉ad(i)，则j-k方向为j→k（否则将存在新的v-结构）。

若DU中包含i→k→j，则i-j的方向为i→j（否则将存在环结构）。

若DU中包含i→k→j和i→I→j，l∉ad(k)，则i-j的方向为i→j（否则将存在新的v-结构或环结构）。

只使用马尔科夫性质，有时不足以得到完整的相依关系，可能存在某些边的方向无法确定的情况，但可作为进一步探索的起点。实际应用中，可以根据数据的具体情况，结合其他方法和经验确定或修改方向。所有边方向确定后，得到有向无环图D对应的链图[3]D*。

步骤3：参数估计。

确定边对应的条件相依关系。结合Vine-Copula结构的参数估计方法，为各边构造藤结构，进行参数估计[3]。本文选取D-vine结构，首先求样本的经验分布函数，用于Copula建模。由于网络结构的方向确定，因此可以确定每个节点的祖先集，对于D中i∈V，j∈ad(i)，k⊆an(i)∩an(j)，根据D-vine的参数估计方法，为边构造一个或多个藤结构，估计所有(xi,xj)|xk的Pair-Copula类型及参数，求得参数对应的Kendall’sτ值，结合实际问题选取|τ|较大的作为该边的相关关系。由于本文主要使用独立性关系作为抽样依据，对于相依关系的依赖较少，因此条件集的选择范围也较宽，在实际问题应用中十分灵活。Pair-Copula函数类型和参数确定后，即可写出相应的联合密度函数，并得到最终的Pair-Copula贝叶斯网络。

2 仿真模拟

本文考虑一个简单的5节点贝叶斯网络，样本容量为5000，使用R语言编写程序。假设给定网络模型以及变量间的相依关系如图1所示。在该网络中，变量间的条件独立关系为2⊥3|1，1⊥4|(2,3)，1⊥5|4，2⊥5|4，3⊥5|4 。首先确定抽样方法，对样本使用上次中的步骤，构建完整Pair-Copula贝叶斯网络。

图1 5节点贝叶斯网络

2.1 抽样方法

根据条件独立关系和相关关系构造变量之间的多项式组，用于抽样。由于任何函数都可用（分段）多项式形式表示，因此可以任意构造能体现相应的独立、相关关系多项式。在多项式回归中，当自变量的幂超过3时，回归系数的解释将变得困难，回归函数也变得很不稳定，对回归模型的应用会受到影响。因而，幂次超过3的模型不常使用。本文以二次方多项式为例，建立如下多项式关系：

抽样步骤如下：

步骤1：根据残差的独立关系，利用独立高斯Copula生成残差ε1,ε2,ε3,ε4。

设ε1,ε2,ε3,ε4为 回 归 残 差 ，ε1⊥ε2，则 (ε1,ε2)～GassianCopula(0,dim=2)，生成两组联合分布为独立高斯Copula的随机数，分别赋予ε1,ε2。生成（0,1）上均匀分布随机数w3～U(0,1)，w4～U(0,1)，根据ε2⊥ε3，有为独立高斯 Copula的条件逆Copula函数。同理

步骤2：根据多项式关系生成样本。首先需生成（0,1）上均匀分布随机数u1。

此时模拟样本得到。

2.2 结构估计

根据1.4中的步骤1，由所有节点的完全无向图出发检验条件独立性，如图2所示。使用基于多项式回归的改进PC算法进行节点历遍后，得到边的具体删留情况如表1所示，此时DAG对应的无向图确定，记为DU。

图2完全无向图

表1 完全无向图边的删留情况

在该例中，5节点结构到此历遍结束。对于更高维的数据，可以不考虑条件集长度大于3的情况，或仅作为参考，综合前几层历遍决定边的删留。

根据1.4中的步骤2，当通过v-结构及其他无向边方向规则为无向图确定方向后，得到链图D*，如图3所示，此时仍有两条边不能完全确定方向，即1—2和1—3。实际问题中可根据情况或结合其他方法确定。暂且假定为1→2，1→3。此时网络完整结构得到，如图4所示。

图3 链图D*

图4有向无环图D

2.3 参数估计

网络结构确定后，并不能最终确定每条边对应的条件相依关系。如边4→5，对应的关系可能是45|12、45|123或其他。根据1.4中的步骤3，为边构造藤结构确定Copula类型，估计参数。本文只考虑Copula函数中最常见的几种函数类型。

最终得到的Pair-Copula类型及参数如表2所示：

表2 Copula类型及参数

联合密度函数为：

Pair-Copula贝叶斯网络得到，如图5所示：

图5 Pair-Copula贝叶斯网络

3 实证

在空调系统设计过程，室外气象参数是负荷计算、设备选型等必需的基础参数，直接影响系统设计和运行效果。气象参数之间具有一定程度的相依关系，本文选取四种典型参数：干球温度、湿球温度、含湿量、太阳辐射，将上文介绍的方法应用到逐时气象数据中进行实证分析，探讨这四种气象参数之间的相关性。数据来自天津市气候中心，选取天津市2010年全年逐时干球温度、湿球温度、含湿量、太阳辐射观测值，每项观测值包含8760个数值。

设干球温度、湿球温度、含湿量、太阳辐射依次为变量x1,x2,x3,x4。根据物理意义及现实情况，温度、湿度等空气状态受太阳辐射强度影响，反之，太阳辐射不受温度等影响，因此确定太阳辐射为首层节点，不考虑在其他变量的条件下太阳辐射与某一变量的关系，即只考虑太阳辐射与其他变量的无条件关系。其他变量的节点顺序未知。

首先使用基于多项式回归的改进PC算法进行节点历遍，选取|τ|=0.3 ，若|τ|＜0.3 认为变量间（条件）不相关，判定满足条件独立性。得到边的具体删留情况如表3所示。

表3 气象数据完全无向图边的删留情况

按照方向规则为无向边确定方向，得到DAG如下页图6所示。

图6气象数据有向无环图

网络结构确定后，将数据变换为Copula水平。由于本文重在研究相关结构，不对边缘分布做具体研究，因此通过经验概率积分变换将逐时干球温度、湿球温度、含湿量、太阳辐射变换为Copula数据。考虑每条边可能对应的所有（条件）关系，估计Pair-Copula类型及参数，并求得参数对应的Kendall’sτ值，结合实际情况和τ值，选取该边的相关关系。最终得到的Pair-Copula类型及参数如表4所示：

表4 气象数据Copula类型及参数

Pair-Copula贝叶斯网络如图7所示。

图7气象数据Pair-Copula贝叶斯网络

从上述Pair-Copula贝叶斯网络看出，干球温度与含湿量、太阳辐射有关；太阳辐射给定时，湿球温度与含湿量相关；含湿量给定时，干湿球温度相关。从物理意义角度分析，干球温度是温度计自由地被暴露在空气中所测量的温度，通常被视作所测量空气的实际温度，与含湿量、太阳辐射有关。湿球温度是温度计的球体表面附着有水时，水份蒸发带走热量后球体的温度。太阳辐射给定时，水的蒸发量跟空气的湿度有关，空气湿度越大蒸发量越小，带走的热量越少，干湿球温度差异越小；空气湿度越小水蒸发量越大，带走的热量也越大，干湿球温差也就越大。通过该Pair-Copula贝叶斯网络可以推断当某一变量发生变化时，干球温度的变化情况。

4 总结

Pair-Copula贝叶斯网络能够有效地转化数据为知识，并利用这些知识进行推理，以解决因果推断等问题。在改进的PC算法的基础上，多项式回归残差的条件独立性检验方法可以良好地检验变量间的条件独立关系，结合Pair-Copula得到完整的相依关系推断模型及联合密度函数。在节点个数较少的情况下，该方法简单有效，适用于在实际应用中解决金融、工程等领域的相关问题。当节点个数较多时，可综合其他方法进一步研究探讨。