一种基于CDC的适用于高维数据的因果推断算法

李洪飞，万亚平,2，阳小华,2，耿家兴

(1.南华大学 计算机学院，湖南 衡阳 421001；2.中核集团高可信计算重点学科实验室，湖南 衡阳 421001)

0 引 言

因果发现旨在开发从观测中学习因果网络结构的算法，而因果结构学习是机器学习和统计领域的新课题之一[1-2]。包括预防科学在内的许多实证科学中，需要研究各种现象背后的因果机制。

通常，从观测数据中发现因果结构是困难的，因为可能的模型的超指数集，造成一些模型可能无法与其他模型区分开来。在变量p数量相当甚至超过样本n数量的情况下是复杂的。尽管该问题本身存在困难，但已经开发了许多算法。目前推断算法一般归为两类：贝叶斯结构学习算法[3-5]和基于加性噪声模型因果发现算法[6-9]。例如PC算法[3]是一种基于约束的贝叶斯结构学习算法，它先推断出一组条件独立关系，然后识别相关的马尔可夫等价类。但这类方法无法区分X→Z→Y和X←Z←Y两种结构。加性噪声模型中，在各种假设下，可从观测数据中恢复准确的图像[10]。如Shimizu等假设数据遵循线性非高斯无环模型(LiNGAM)[7]，后续的DirectLiNGAM[8]和PLiNGAM[11]方法通过两两统计量迭代选择因果顺序。LiNGAM提出的这些方法在高维设置中是不一致的，这些设置允许变量数的伸缩与样本大小相同或更快，并且只在线性模型下适用。Hoyer针对非线性数据因果模型，提出适用于连续数据算法ANM[6]。为扩展到离散数据的情况，对ANM算法进一步改进，提出后非线性数据的PostNonLinear算法[9]和信息几何方法[12]。

近年来，Makoto等提出最小二乘独立回归算法(LSIR)[13]。该算法通过交叉验证自然地优化内核宽度和正则化参数等调优参数，从而避免以依赖数据的方式过度拟合。这些方法都只考虑了少量的变量，一旦在高维的情况下(n>7)，因果发现方法的输出可能高度依赖于顺序，准确率会很低。在大多数情况下是两个变量进行方向识别。因为根据理论和算法，可扩展性有限，所以这些因果推断方法都不适用在高维数据的网络结构中。

为解决在高维数据下学习到较准确的因果网络结构问题，曾千千等提出基于MIC和MI的因果推断算法[14-15]。这些算法将高维网络结构学习问题分解成网络中每个节点对应的低维因果网络结构学习问题。这些算法中，基于信息论的互信息(MI)或者最大信息系数(MIC)能够较好地删除不相关的节点，降低计算复杂度，但MIC对于高维变量不可用，对于大样本问题耗时MI也不稳定。随着样本量的增大，方差没有变小，对异常值的鲁棒性较弱。由Jiang等提出的CDC相比较互信息和最大信息系数，可更为准确地检测出变量之间的关联关系[16]。

基于此，文中提出一种基于CDC的高维数据因果推断算法。该算法利用CDC对变量间的依赖度进行检测，再利用条件独立性检测精炼目标节点的父子节点集，然后使用非线性最小二乘独立回归(LSIR)(两个变量之间效应方向的新测度)，为图中的目标节点及其父子节点之间的无向边标注因果方向，最后迭代所有的节点完成完整的因果网络结构。由于CDC对大样本下异常值的鲁棒性强，因此提高了算法的准确率。实验也表明在高维数据下该算法的精确度优于其他因果推断算法。

1 预备知识

1.1 加性噪声模型

1.2 最小二乘独立回归

最小二乘独立回归(LSIR)是通过最小化输入和残差之间的平方损失互信息估计值来学习加性噪声模型，识别两个非线性变量的因果关系的方法。LSIR相对现有方法的一个显著优势是，调优参数(如内核宽度和正则化参数)可以通过交叉验证自然地进行优化，从而避免以依赖数据的方式过度拟合。LSIR与最先进的因果推理方法相比是比较有利的。

(1)

采用平方损失相互信息(SMI)[14]作为独立措施:

(2)

使用SMI，获得一个分析形式的估计值。得到SMI估计量后，开始学习参数β回归模型：

这种方法称为最小二乘独立回归。

对于回归参数学习，可以简单地采用梯度下降法:

(4)

其中，η是步长，可以选一些近似线搜索方法，如Armijo规则[14]。

(5)

(6)

为了确定因果关系的方向，计算两个方向X→Y(即X引起Y)和X←Y(即Y导致X)的p值pX→Y和pX←Y。对于一个给定的显著性水平δ，确定的因果方向如下。

如果pX→Y>δ和pX←Y≤δ，模型选择X→Y；

如果pX←Y>δ和pX←Y≤δ，X←Y模型被选中；

如果pX→Y,pX←Y≤δ，X和Y之间没有因果关系；

如果pX→Y,pX←Y>δ,则建模假设是不正确的。

总之，检验因果模型Y=fY(X)+EY或替代X=fX(Y)+EX是否与数据吻合较好，其中拟合优度是通过输入与残差之间的独立性来衡量的(即噪声估计)，输入和残差的独立性可以通过置换测试在实践中确定[14]。

1.3 Copula相关系数

Copula相关系数(CDC)用来测量两个随机向量X和Y之间的依赖关系，是一种较好的关联检测鲁棒依赖测度。CDC对异常值更具有鲁棒性，适用于更广泛的模型，尤其适用于高维问题。

定理1(概率积分变换)：对于具有累积分布函数(CDF)F的连续随机变量X，随机变量U=F(X)在[0,1]上均匀分布。因此，向量：

U=(U1,U2,…,Up)=P(X)=(F1(X1),F2(X2),…,Fp(Xp))

(7)

称为Copula变换，具有均匀的边缘。

定理2：让随机向量X=(X1,X2,…,Xp)连续边际CDFs,Fi,1≤i≤p。然后X的联合累积分布函数F(X),X=(x1,x2,…,xp)，唯一表示为:

F(X)=Cx(F1(x1),F2(x2),…,Fk(xp))

(8)

其中分布函数Cx称为X的Copula。

定理3：设X和Y为连续随机变量，那么X和Y是独立的，当且仅当CXY(FX(x),FY(y))=FX(x)FY(y)，其中FX(x)和FY(y)分别为x和y的分布函数。

首先考虑特殊情况p=q=1。如果X和Y是独立的，将FX(X)andFY(Y)分别表示为X和Y的分布函数，由定理3可知：

CXY(FX(x),FY(y))=FX(x)FY(y)

(9)

如果X和Y是独立的，有W=FX(X)andV=FY(Y)是独立的。因此，将多维依赖检验问题转化为零假设(H0)，W和V是独立的二维依赖检验问题，可以通过MIC、互信息等提出的另一种依赖测度来求解。

定义1：设X和Y为两个随机变量，X和Y的最大相关系数(MCC)为：

(10)

如果将φ1,φ2限制为线性函数，MCC是皮尔森相关系数。为了计算MCC，通常会将限制φ1,φ2限制为再生核希尔伯特空间理论(RKHS)，例如n阶的多项式函数的空间。而ACE用来计算MCC。

定义2：利用上述表示法，两个随机向量X和Y之间的CDC由FX(X)andFY(Y)之间的MCC给出：

CDC(X,Y)=MCC(FX(X),FY(Y))

(11)

其中，FX(x)andFY(y)分别为X和Y的分布函数。

可以看出，CDC是一个基于秩的依赖性测度，因为分布函数是基于秩的，所以CDC适用于高维变量。在实际应用中，真实的边际分布函数是未知的，用经验边际分布或估计边际分布代替，得到CDC的估计。

2 因果结构学习算法

本节根据相较于互信息和最大信息系数更为准确地检测出变量间关联关系的CDC，结合发现因果骨架的结构学习，二元变量的方向识别，最终迭代得到一个完整表现出高维数据间因果关系的因果网络结构。研究表明，张等基于互信息与曾等基于最大信息系数的算法构成的无向图，对大样本和高维下异常值的鲁棒性不高，一定程度会增加父子节点集，造成条件独立测试的计算复杂度，得到的最终网络与真实网络差异很大[14-15]。相较于MIC和MI，CDC能够更好地度量节点间的依赖关系，减少冗余带来的计算复杂度，从而提高效率，并获得更加准确的因果网络结构。该算法基本框架如图1所示。

图1 算法的基本框架

2.1 无向子图构造

CDC可以快速对随机变量之间的依赖关系进行检测。在样本量增多的情况下，CDC是最佳的依赖测度，它对异常值具有鲁棒性，计算效率高，对大多数函数类型具有较好的适用性。对随机变量X和Y计算CDC(X,Y)，X和Y之间的依赖程度与CDC的值成正比。当CDC值越大，则两个变量对应的点相连程度越高。反之，若X和Y之间CDC值很小，甚至达到0，则两个变量之间是独立的，对应的两点不相连。根据CDC的优点计算变量间的CDC(X,Y)值构造无向图，步骤如下：

(1)对数据集中任意的两个随机变量计算其之间的CDC值，找到每个节点与其他节点的CDC值中最大的值，用MCDC(X,Y)表示；

(2)对于任意节点，假设为y，其他节点表示为X集，如果满足下列条件，则直接连接在变量X和Y之间建立：CDC(xi,y)≥α·MCDC(X)，否则将xi从X中移除；

(3)将任意的xi∈X加入到PC(y)，应用条件独立性测试(CI)删除PC(y)中非父子节点；

(4)重复步骤3，迭代完X中所有节点。

上述步骤中，0.5≤α≤1是一个参数，来确定此阶段的连接数量，如果该参数非常高(接近1.0)，在该算法的这个阶段，存在多个真边被拒绝的高概率。另一方面，将该参数设置的低，可能会导致在算法的这个阶段包含几个错误的边。对于该研究中的问题，参数α设置为0.9，并且发现包括大多数真实边缘，同时允许在网络结构中仅添加少量错误边缘。

2.2 目标节点和父子节点集PC(y)中每个节点进行方向识别

在前一阶段用CDC、CI测试得到的父子节点集，构成了一个准确的骨架。再利用LSIR算法对骨架之间每两个节点之间的无向边进行方向识别。相比较于ANM、IGCI二元算法，LSIR在大样本下通过交叉验证进行优化，避免了过分依赖数据而导致的过度拟合。虽然LSIR算法不能有效处理高维数据下的因果网络结构学习问题，但由于在算法2.1无向图构造过程中，已经得到了目标节点y的父子节点集合PC(y)，所以可以利用LSIR对PC(y)中的每一个变量和y之间进行方向判别，这等同于在2维间应用LSIR。这种分治策略保证了它的有效性。

2.3 算法描述

文中算法记为CDC & LSIR causal inference (CLCI)。

Input:variable set X={x1,x2,…,xn};threshold αOutput:DAG G1 for each xi∈X doSet xi=y,xi∈X,PC(y)={}/∗ 选取X中任意一个节点为目标节点∗/2 for each xi∈X do /∗ 对PC(y)进行精炼简化,除掉非父节点∗/if CDC(y;xi)<α∗MCDC(X), then X=Xxi3 for each xi∈X do /∗ 构造关于节点y的无向图∗/ add xi to PC(y),if there is a set S ( arbitrary subset of PC(y)xi),such that y⊥xi|S,then PC(y)=PC(y)xi4 for each xi∈PC(y) doif there is a set S(arbitrary subset of PC(y)xi), such that y⊥xi|S, then PC(y)=PC(y)xi5 for each xi∈PC(y) do /∗对骨架之间每两个点的方向进行判别∗ /employs LSIR algorithm to distinguish the direction of (y,xi)

3 模拟数据实验结果与分析

3.1 实验设置与环境

通过人工合成模拟数据进行实验，实验一用于在大样本下测试CDC的性能，而实验二用于验证CLCI算法的有效性。实验中尽可能假设复杂条件，即观察样本数目足够大，并且模糊了对噪声变量非高斯性的假设。每个节点在因果网络结构中的序列按照函数y=ω1*tan(x1)+ω2*tan(x2)+ε生成。每个函数中的权值分别为ω1、ω2,随机取值0.3～0.7之间，y的父节点分别为X1和X2。噪声ε服从均匀分布，且权值为6%。文中提出的算法记为CLCI。实验平台为Window 7 64 bit操作系统，配置为Intel酷睿i5-4200U，内存4 GB，实验环境为Matlab-R2016a。

3.2 实验分析

在不同样本下，CDC与MI和MIC计算时间对比如图2所示(表示时间效率)。

图2 MI、MIC和CDC在不同大样本下的计算时间对比

从结果可以看出，随着样本量的增大，MIC的运行时间迅速增加，MI直到样本大小约为4 500时才开始缓慢增加，而CDC的运行时间则没有增加趋势。

为评价CLCI在高维下的优劣，避免随机性，在实验中选择5 000样本量分别生成20维到200维的因果网络图，用来测试高维度下不同样本量的实验效果。实验引入三个常用指标recall、precision和F1来评价算法性能。在不同维度下四种算法的评分参数如图3所示。

从图3可以看到，IGCI算法的召回率在超过40维左右时比其他三种算法都要高，但由于在高维因果网络结构下，IGCI错误地添加了很多边，准确率相对较低。同样，加性噪声模型(ANM)算法在高维数据情况下三个评分参数比CLCI和TPCI算法低。由于CLCI和TPCI算法采用了降维的方法，在维数增加的过程中，其保持了较好的稳定性，但是相对于利用MI的TPCI，CLCI利用CDC在大样本下的鲁棒性更好，时间复杂度更低，耗时更短。由图3也可以看出，CLCI要优于其他三种算法。

(a)recall

(b)precision

(c)F1

4 结束语

与其他适用于高维数据的因果推断方法不同，文中结合基于信息论的Copula依赖系数，可以在高维和大样本数据下更为准确地检测出变量之间的关联关系，降低计算复杂度，再结合条件独立性测试对数据集进行无向结构学习，最后用LSIR算法对无向结构骨架进行每两点间的方向判断，得到最终的因果网络结构。文中采用虚拟网络进行的实验表明，利用CDC进行无向骨架的构造，耗时短，计算复杂度低于其他算法。当数据集维数较高、样本量大时，利用分治策略，该算法要大大优于其他因果推断算法。