郝志峰,丁凯培,蔡瑞初＊,陈薇
(1. 广东工业大学计算机学院,广东 广州 510006;2. 汕头大学理学院,广东 汕头 515063)
挖掘可观测的时间序列之间的因果关系可以帮助人们理解复杂事物的本质和规律。近年来,因果关系被广泛应用于深度学习[1]、生物信息学[2]、金融经济[3]、社会科学[4-5]等领域。通常来说,随机实验是挖掘因果关系的黄金标准[6],但受限于经济成本和道德伦理方面的约束,随机实验并不总是可行的。因此,如何从可观测的时间序列中挖掘因果关系成为数据科学的主要问题。
当前,研究人员已经提出了多种恢复时序数据集的因果结构的方法,常见的有基于约束的方法[7-9]、基于因果函数模型的方法[10-11]等。然而,大多数方法都未考虑数据集中非稳态扰动带来的影响,进而导致方法失效。非稳态扰动引起的数据分布偏移常见于真实的应用环境中,由非稳态扰动所导致的统计误差会一定程度上打乱数据分布,使数据分布难以表征。一些学者尝试用数学工具刻画非稳态扰动的规律,如文献[12]提出基于约束的异构/非平稳因果发现(CD-NOD)算法,该算法通过引入代理变量从而对非稳态扰动进行表征,并将包含代理变量的增广变量集作为因果发现算法的输入。CD-NOD算法首先利用PC(Peter-Clark)算法检测局部因果机制发生变化的变量,并恢复变量的因果骨架图,然后通过多项定向规则对因果骨架图中的无向边进行定向。然而,CD-NOD算法额外假设数据集中非稳态扰动是彼此独立的,该假设致使CD-NOD算法无法从不满足非稳态扰动独立的数据集中恢复因果结构,并且存在马尔可夫等价类问题。
针对以上问题,本文基于非稳态扰动与时序信息存在高度相关性的规律,在加性噪声模型的基础上将非稳态扰动刻画成一项关于时序信息的映射并融入模型中,从而提出一种非稳态加性噪声模型(Nons-ANM)。给出该函数模型的识别条件,并基于模型的识别条件提出一种两阶段的因果发现算法。第1阶段基于因果结构图中叶子节点的噪声与其余节点均独立的性质找到因果次序;基于前阶段的因果次序先验,第2阶段进一步移除冗余的父子关系,最终得到变量集的因果结构。
针对观察数据特性的不同,因果发现方法又可以分为基于时序观察数据的方法和基于非时序观察数据的方法[13]。
在非时序方法中,基于约束的方法通常利用条件独立和忠实性假设来恢复变量之间的因果结构,经典的方法包括文献[14]介绍的PC算法和文献[15]介绍的IC(Inductive Causation)算法。这类方法无法解决马尔可夫等价类的问题,因此算法最终得到的因果结构图中可能包含未确定方向的边。为解决马尔可夫等价类的问题,学者们提出了基于函数的方法。这种方法的核心思想在于:假设存在一种表达力较强的函数形式能够刻画变量之间的关系。在这种函数关系下,结果变量能够由原因变量和噪声变量的某种映射所表征。在一定的约束条件下,变量之间的因果关系能够通过检验原因变量和噪声变量的独立性得出。经典的基于函数的方法包括文献[16]介绍的非线性加性噪声模型、文献[17]介绍的线性非高斯无环模型和文献[18]介绍的后非线性因果模型等。
时序观察数据在时间维度上伴随着“原因在前,结果在后”的特性[13],基于这项特性,非时序方法能够向时序数据的场景拓展。典型的有文献[7]提出的PC瞬时条件独立性(PCMCI)算法,该算法在PC算法的基础上拓展出恢复多元时序数据集的因果结构的能力。具体来说,PCMCI算法首先使用 PC算法学习得到数据集的马尔可夫等价类集合,然后再利用瞬时条件独立性(MCI)检验以降低因果关系的误发现率。然而,PCMCI算法存在未考虑即时效应和隐变量的局限性。因此,文献[8]提出PCMCI+算法以解决即时效应的问题;文献[9]提出的LPCMCI算法进一步解决了隐变量的问题。类似地,文献[10]提出的TiMINo算法和文献[11]提出的VarLiNGAM算法则是基于函数模型的方法在时序场景上的拓展。
当前针对非稳态数据的因果发现算法的研究仍处于初始阶段,文献[12]介绍的CD-NOD框架着力于在变量的分布发生变化的情况下恢复因果结构,该框架的关键部分为引入代理变量C充当潜在变化因素的表征。CD-NOD能够恢复非稳态数据的因果骨架,但是非稳态扰动的独立变化假设和马尔可夫等价类难题限制了该算法对完整因果结构的学习能力。具体来说,马尔可夫等价类是指具有相同的因果骨架和条件概率分布,但无法确定方向的结构。在CD-NOD算法的定向规则失效后,该算法最终仅能输出由PC算法学习到的部分定向的因果结构。
本文所考虑的问题为:给定非稳态时间序列数据集D,如何有效地根据非稳态扰动的特征建模,并恢复非稳态数据集D中n个变量的因果结构。
针对所研究的问题,本节提出一种非稳态加性噪声模型,并给出该模型在非稳态场景下的可识别性条件及其证明,随后基于非稳态加性噪声模型的可识别性理论提出一种两阶段的因果关系发现算法。
3.1 非稳态加性噪声模型
在稳态场景中,文献[16]基于二元变量的因果关系提出非线性加性噪声模型,并证明该模型在限定的函数和噪声分布下具有可识别性。文献[10-19]则进一步将非线性加性噪声模型推广至时序数据场景和多元变量场景。然而,目前对于非稳态加性噪声模型的建模及其识别性仍然未知,本文首先给出非稳态加性噪声模型的定义,随后给出模型的识别性条件及证明。
定义1(非稳态加性噪声模型) 给定一组观测变量集X={X1,X2,…,Xn},若这组观测变量满足以下约束,则称为非稳态加性噪声模型:
1)观察变量具有一定的因果次序,因果次序在后的变量不能是次序在前的变量的原因;数据产生过程是递归的,变量间的因果关系可以表示为一个有向无环图。
2)伪因果充分性[12]假设:若存在未观测到的混淆变量,则这个混淆变量可以用关于时序信息T的光滑函数表示。在任意时刻T=t,这项混淆变量的函数值是确定的。
3)观测变量Xi∈X服从以下函数关系:
(1)
受文献[19]提出的多元加性噪声模型的可识别性条件启发,Nons-ANM的可识别性理论如定理1所示。在给出定理1之前,先引入(B,F)-IFMOC的识别理论。
引理1(B,F)-IFMOC的识别理论[19-20]:假设一组变量X的联合分布p(X)是由以因果图G表征的(B,F)-IFMOC的函数模型导出的,那么这组变量一定不存在两个不同的(B,F)-IFMOC函数模型表征,同时也一定不能由两个不同的因果图G和G′表征。
定理1假设变量集X={X1,X2,…,Xn}满足非稳态加性噪声模型,模型属于(B,F)-IFMOC空间,且增广变量集X∪{T}的因果结构图服从因果最小性假设,那么变量集X的因果结构能够被唯一确定。
证明假设增广变量集X∪{T}存在两种不同的非稳态加性噪声模型表征,分别对应两个不同的因果结构图G和G0。
定理1提供了非稳态加性噪声模型的识别性条件,而对于模型中存在直接因果关系的变量对之间的非对称性质,可由定理1推论而来。
证明在本文模型中,若两个变量之间存在因果关系,如Xi→Xj,可分为两种情况:1)存在时滞的因果关系;2)即时因果关系。两种情况不会同时发生,因此分别证明在以上两种情况下本文模型的非对称性质均成立。
推论1提供了本文模型中原因变量和结果变量之间的非对称性质。基于上述理论,第3.2节将提出面向非稳态数据的因果发现算法。
3.2 非稳态数据的因果发现算法
基于定理1提供的非稳态加性噪声模型的识别性条件,本节提出一种两阶段的因果发现算法。算法的初始输入为观测变量的时序数据集,最终输出为表征观测变量因果结构的有向无环图,算法框架如图1所示。
图1 Nons-ANM算法框架Fig.1 Nons-ANM algorithm framework
算法第1阶段主要基于有向无环图的汇点的性质:对于一个有向无环图,必然存在出度为0的叶子节点,本文将其称为“汇点”。在因果结构图中,若Xsink为汇点,则其噪声变量为Nsink,它将拥有以下性质:
Nsink⊥{Xk|Xk∈X{Xsink}}
(2)
对于因果结构图中的汇点Xsink,其噪声Nsink将独立于变量集中其余所有变量,这是因为其余的变量Xk均为汇点Xsink的非后裔节点,因此可由不包含汇点Xsink成分的独立同分布噪声变量集N
基于上述性质,变量的因果次序可通过迭代寻找汇点的方法来获得。具体来说,对于未确定因果次序的变量Xs∈X,利用变量集中其余未确定因果次序的变量{Xk|Xk∈X{Xs}}对Xs进行回归并计算得到残差rs,随后计算该残差与其余变量的独立性得分IND(rs,{Xk|Xk∈X{Xs}}),该得分衡量了变量间的独立性强弱,衡量方法IND(·)可根据数据性质灵活选择,本节列举了一些参考方法。重复上述步骤,选出独立性得分最高的变量作为汇点,将汇点插入因果次序中,并将汇点的数据从数据集中移除。算法伪代码如算法1所示。
算法1因果次序学习
输入观测变量的时序数据集D={x1,x2,…,xn},最大时延p
输出观测变量的因果次序Π
1.初始化变量索引集S={1,2,…,n},因果次序Π为空集
2.WHILE !is Empty(S):
3.FOR s in S:
5.计算IND(rs,{xk|k∈S{s})并记录
6.END FOR
7.选出独立性得分最高的变量XM,从变量索引集S中移除M,在因果次序Π中添加M
8.END WHILE
9.RETURN Π
算法1可以学习得到变量集X的因果次序Π,基于因果次序Π可以生成变量集X的部分因果结构。具体来说,在因果次序中,次序在后的变量不能是次序在前的变量的原因。这就意味着,因果次序中包含的确定信息为不存在由次序在后的节点指向次序在前的节点的有向边;而未确定的信息为次序在前的节点可能存在指向次序在后的节点的有向边,同时也可能不存在有向边。因此,因果次序仅表征了因果结构的部分信息,需要进一步消除冗余信息,从而得到最终的因果结构图。
消除冗余因果关系的步骤如下:
1)对于各节点Xi∈X,以Xi为结果变量,以因果次序在Xi之前的全部节点为原因变量,连接原因变量和结果变量,初始化变量集X的部分因果结构G*。
2)对于G*中可能存在的冗余边,采用回归计算得到残差然后进行独立性检验的方式消除,即对于各节点Xi∈X,依次假设因果次序在Xi之前的某个节点Xj∈{Xj|Xj∈X{Xi},Π(j)<Π(i)}不是Xi的直接原因;在特征数据集中移除Xj的数据后再对Xi的回归计算得到残差ri,随后检验残差ri与变量集{Xj|Xj∈X{Xi},Π(j)<Π(i)}的独立性。若独立性仍然成立,则认为Xj不是Xi的直接原因,并移除冗余边Xj→Xi;否则认为Xj是Xi的直接原因,并保留连接边。算法伪代码如算法2所示。
算法2移除冗余因果边
输入因果次序Π,观测变量的时序数据集D={x1,x2,…,xn},最大时延p,显著性水平α
输出因果结构图G
1.基于因果次序Π初始化有向无环图G*。
2.FOR i in Π:
3.FOR j in {j|Π(j)<Π(i)}:
5.IF IND(ri,{xk|Π(k)<Π(i))>α:
6.删除G*中的有向边Xj→Xi,更新父亲节点集PAi=PAi{j}
7.END IF
8.END FOR
9.END FOR
10.RETURN G*
根据数据集特征,算法中的回归方法和独立性检验方法可以相应调整。在因果发现中常见的回归方法有高斯过程回归(GPR)、向量自回归模型(VAR)等;常见的独立性检验方法有偏相关检验[22]、希尔伯特-施密特独立性准则(HSIC)[23-24]、互信息(Ml)[25]、核条件独立性检验(KCI-test)[26]等。本文算法的实际表现受回归工具和独立性检验工具的影响,统计学的难题会影响算法的实际表现,例如高维变量集、小样本量等场景。
对本文算法进行复杂度分析:算法1和算法2的时间复杂度均为O(n2·Or(n,l)·Oi(n,l)),其中,n为观测变量的个数,l为时序数据集的样本量,Or(n,l)表示算法选用的回归方法的时间复杂度,Oi(n,l)表示算法选用的独立性检验方法的时间复杂度。
本节将分别基于仿真数据集和真实因果结构数据集开展随机实验,从而对本文算法的有效性进行评估。另外采用CD-NOD算法、LPCMCI算法和TiMINo算法作为对比方法,其中CD-NOD算法为非稳态因果关系学习算法,LPCMCI算法是学习含有即时因果效应和隐变量的因果关系学习算法,TiMINo算法则是经典的基于函数因果模型的时序因果关系学习算法。
仿真数据集实验和真实因果结构数据集实验中的对比算法均使用官方实现的开源代码,算法参数均设置为默认值。本文算法Nons-ANM采用高斯过程回归作为回归方法,采用基于希尔伯特-施密特独立性准则的假设检验作为独立性检验方法。实验中所有方法的显著性水平α均设置为0.05。
在实验部署环境方面,仿真数据集实验和真实因果结构数据集实验均执行于内存为8 GB、CPU型号为Intel i7-7700、操作系统为Windows的环境。
4.1 仿真数据集实验
在仿真数据集实验中,本文基于Erdos-Renyi模型随机生成不同的有向无环图来表示因果结构,并设计了5组控制变量实验,具体设置为:变量数量为{5,10, 15, 20},平均入度为{1.0,1.5, 2.0, 2.5},样本数量为{1 000,2000, 3 000, 4 000},非稳态变量的占比为{0.0,0.4, 0.8, 1.0},最大时延为{0, 2, 4},花括号中的加粗参数为默认的参数设置。变量集中的变量分为稳态和非稳态两类,对于稳态的变量,其时序数据基于函数因果模型生成,如式(3)所示:
(3)
对于非稳态扰动的变量,其时序数据的生成则基于函数因果模型,如式(4)所示:
(4)
实验的评价指标选用准确率(P)、召回率(R)和F1值(F1),各评价指标的计算公式如下:
(5)
(6)
(7)
其中:TP表示学习正确的连接边的数量;FN表示原结构不存在边,且算法结果中也没有边的数量;FP表示原结构不存在边,但算法结果存在边的数量。在这套评价算法中,正确率是指学习的因果结构图中正确的边数占比;召回率表示正确识别的边数在所有真实存在的正确边数中的占比;而F1值结合正确率和召回率,能够综合地衡量一种算法的有效性。
本文算法与CD-NOD、LPCMCI和TiMINo算法的实验结果如图2所示。
图2 各算法在不同参数下的实验结果Fig.2 Experimental results of each algorithms in different parameters
1)图2(a)～图2(c)给出不同变量数量下的实验结果。可以看出,本文算法在所有变量数量的设定下均取得最高的准确率、召回率和F1值,这验证了本文算法的鲁棒性。CD-NOD算法的实验结果呈现出准确率高于召回率的特点,原因是CD-NOD算法的定向规则中的独立扰动假设造成了定向误差。LPCMCI算法拥有从包含隐变量的数据集中恢复因果结构的能力,这项能力帮助LPCMCI算法在实验默认设置下取得了0.56的F1值,但其总体实验结果略低于CD-NOD算法,这是因为LPCMCI算法针对隐变量的能力并不完全适用于非稳态扰动的场景。TiMINo算法采用了保守的策略,当算法无法确定因果关系时将执行熔断机制;由图3(a)可知,TiMINo算法在实验默认设置下熔断率为0.8;而当变量数量增至20时,其熔断率达到1.0,即TiMINo算法无法恢复任何一组仿真数据集的因果结构。由此可知,经典的时序因果发现方法难以直接应用于非稳态场景。
图3 TiMINo算法在不同参数设置下的未识别率Fig.3 Unrecognition rate of TiMINo algorithm under different parameter settings
2)从图2(d)～图2(f)可以看出,实验中的所有算法对样本数量都比较敏感。在样本数量逐渐增加的过程中,本文算法的F1值逐渐提升,这是因为高样本数量能降低统计误差,从而帮助本文算法逼近理论效果。在所有样本数量的实验设置下,本文算法均取得了最高的F1值,这体现了本文算法的优越性。
3)图2(g)～图2(i)展示了不同平均入度下的实验结果。可以看出,CD-NOD算法和LPCMCI算法对平均入度较敏感,两个算法的F1值随着平均入度增加而明显下降,这表明稠密的因果结构更易于暴露算法的不足。TiMINo算法的熔断率在所有设定中均不小于0.7,这表明了TiMINo算法对非稳态数据集的学习能力不足。相比之下,本文算法在所有设置下的均表现稳定且优势明显,这表示本文算法能够有效地学习稠密的因果结构。
4)图2(j)～图2(l)给出了不同非稳态节点占比的实验结果。可以看出,在所有算法中本文算法拥有最好的非稳态因果关系学习能力,即使在非稳态节点占比为1.0时,本文算法依然能取得0.81的F1值。随着非稳态扰动节点占比的增加,其余3个算法的实验结果均大幅降低,这进一步体现了本文算法的鲁棒性和优越性。
5)图2(m)～图2(o)给出了不同最大时延下的实验结果。最大时延的增加会增大时序因果网络的复杂性,进而增大因果发现算法的学习难度。从图中可以看出,CD-NOD算法、LPCMCI算法和TiMINo算法在最大时延为4时难以取得有效的实验结果,这是因为以上3种算法均没有学习复杂的时序因果网络的能力。而本文算法在所有最大时延的实验设置下表现稳定,实验结果显著优于其他3种算法,说明了本文算法的有效性。
6)图3表示TiMINo算法在不同控制实验中触发熔断策略的数据集比率,即TiMINo算法的未识别率。可以看出,TiMINo算法频繁触发了熔断机制:在所有控制实验中,TiMINo的未识别率均不小于0.7;而在部分复杂的因果结构中的未识别率达到1.0,例如20个节点和非稳态节点占比达到1.0等。这进一步表明,未考虑非稳态扰动的时序因果发现算法难以直接应用于非稳态场景中。
4.2 真实因果结构数据集实验
本节采用3种不同规格的真实因果结构进行实验,相关的真实因果结构信息可参考开源网站:https:∥www.bnlearn.com/bnrepository/。本节实验中选用的真实因果网络的参数表1如所示。
表1 真实因果网络结构参数Table 1 Real causal network structure parameters
本节实验中采用的算法参数和实验参数均与仿真数据集实验相同,同样采用准确率、召回率以及F1值作为算法效果的评价指标。
各算法的实验结果如图4所示。由图4可知,本文算法在所有数据集中均取得了最优的实验结果。特别地,本文算法在sachs结构的数据集中取得的F1值为0.86,相比CD-NOD的0.52提升了65.4%,相比LPCMCI的0.53提升了62.3%,这是因为sachs结构比较稠密,而本文算法在稠密的因果结构下更加能够体现其优越性。而对于变量个数少且相对稀疏的因果结构(如survey结构),本文算法和其他对比算法的F1值接近,这则是因为简单的因果结构难以体现非稳态扰动造成的影响。总体而言,真实因果结构数据集的实验结果验证了本文算法的有效性。TiMINo算法在不同真实结构下的未识别率如图5所示。
图4 不同算法的真实结构实验结果Fig.4 Experimental results of real structures of different algorithms
图5 TiMINo算法在不同真实结构下的未识别率Fig.5 Unrecognition rate of TiMINo algorithm under different real structures
现有因果发现方法难以恢复非稳态数据集的因果结构,本文提出非稳态加性噪声模型,并基于该模型的识别性理论提出一种鲁棒的两阶段算法,从而为非稳态数据集的因果结构学习难题提供解决方案。为验证本文算法有效性,引入CD-NOD算法、LPCMCI算法和TiMINo算法,分别在仿真数据集和真实结构数据集上部署实验,最终的实验结果体现了本文算法相较于主流因果发现算法在准确率和鲁棒性等方面的优势。下一步将对模型的加性噪声假设进行弱化并优化算法流程,以扩展本文方法的适用范围。
猜你喜欢 加性次序时序 《汉纪》对汉帝功业次序的重构及其意义历史教学问题(2023年6期)2024-01-18基于时序Sentinel-2数据的马铃薯遥感识别研究中国农业信息(2023年3期)2023-03-18ℤ2ℤ4[u]-加性循环码阜阳师范大学学报(自然科学版)(2022年1期)2022-04-02基于Sentinel-2时序NDVI的麦冬识别研究中国农业信息(2021年3期)2021-11-22企业家多重政治联系与企业绩效关系：超可加性、次可加性或不可加性系统管理学报(2018年3期)2018-08-13企业家多重政治联系与企业绩效关系：超可加性、次可加性或不可加性系统管理学报(2018年2期)2018-08-13生日谜题小天使·六年级语数英综合(2017年8期)2017-08-04一种毫米波放大器时序直流电源的设计电子制作(2016年15期)2017-01-15基于加性指标的网络断层扫描的研究哈尔滨师范大学自然科学学报(2015年1期)2015-04-19放假一年小天使·五年级语数英综合(2014年11期)2014-11-06