因果算法伪代码

一、LiNGAM线性非高斯模型的伪代码

1. DirectLiNGAM 直接线性非高斯模型

DirectLiNGAM是直接线性非高斯模型，其第一步是基于互信息对特征进行排序，第二步是基于特征排序结果获得特征间的因果关系强度。DirectLiNGAM对m*n的矩阵，时间复杂度O(m*n^3)，1000*1000矩阵要跑1000×170秒，速度很慢。

伪代码：

特征排序：

这一步的目的是获得特征的排序结果K，例如K=[x1,x2,..,xn],那么因果的方向只能从左到右，不能从右到左，这样就是一个有向无环图的框架。x1-->x3是合法的，x3-->x1是非法的，xi的原因只能是x1,x2,...,x(i-1)的线性组合

X为m*n的数据，m是样本数，n是特征数。

U=所有特征

K=[]

X_=copy(X)

while U!=[]:

1.1 计算每个特征的互信息mutual_info，把互信息mutual_info最小的特征记为m，过程如下：

for i in U:

mutual_info(i)=0

for j in U:

if i!=j:

std(x)=(x-x.mean())/x.std()

mutual_info(i,j)=entroy(std(X_[:,i]))-entroy(std(X_[:,j]))

mutual_info(i)+=min(0,mutual_info(i,j))**2

1.2 对特征m根据残差对偶函数更新矩阵：

for i in U:

if i!=m:

残差对偶函数：residual(xi,xj)=xi-(cov(xi,xj)/var(vj))*xj

根据残差对偶函数更新矩阵： X_[:,i]=residual(X_[:,i] ,X_[:,m])

1.3 m加入K中并从U中删除m：

K.append(m)

U.pop (m)

计算因果权重：

记B是n×n的矩阵，初始化为0， Bij 表示xj到xi的因果关系强度

for i in range(1,len(K)):

2.1使用K[:i]（K 的前i-1个特征）预测K[i]（第i个特征），也就是构造X[:,K[:i]] （n×（i-1））到 X[:,K[i]] （n×1）的映射，并返回参数coef，方式如下：

用LinearRegression拟合 X[:,K[:i]] 到X[:,K[i]]的映射，并返回特征的权重系数weight_lr

用 LassoLarsIC拟合 X[:,K[:i]]×weight_lr 到X[:,K[i]]的函数映射，并返回参数weight_lasso

最终coef=weight_lr×weight_lasso

2.2更新因果关系到矩阵B

B[K[i],K[:i]]=coef

2. ICALiNGAM 线性非高斯模型

ICALiNGAM 、DirectLiNGAM都是线性非高斯模型，但ICALiNGAM是基于ICA矩阵分解获得混淆矩阵对特征进行排序的。

伪代码：

特征排序：

这一步的目的是获得特征的排序结果K，例如K=[x1,x2,..,xn],那么因果的方向只能从左到右，不能从右到左，这样就是一个有向无环图的框架。x1-->x3是合法的，x3-->x1是非法的，xi的原因只能是x1,x2,...,x(i-1)的线性组合。

X为m*n的数据，m是样本数，n是特征数。使用ICA做特征排序的过程如下：

1.1 应用ICA算法获得一个矩阵分解 X=SA ,A是一个特征间的混淆矩阵（A的大小为n*n），我们只处理W= A^(−1)。

1.2 对W做行变换产生一个矩阵 W1 ，其对角线的元素都是非零的。由于微小的估计误差将导致W的所有元素非零, 因此行 W1 是满足最小化 Σi 1/ |W1(ii)| 的解

1.3 将W1的每行除以其相应对角线上的元素，获得一个新的对角线元素全为1的矩阵W2

1.4 获得估计矩阵C = I - W2

1.5 找到C的置换矩阵P (对行、列做相同的初等变换)，产生一个对应的矩阵

C1 = PC(P^T),使得C1是严格的下三角矩阵，从而根据P获得相应的因果顺序K

计算因果权重：

记B是n×n的矩阵，初始化为0， Bij 表示xj到xi的因果关系强度

for i in range(1,len(K)):

2.1使用K[:i]（K 的前i-1个特征）预测K[i]（第i个特征），也就是构造X[:,K[:i]] （n×（i-1））到 X[:,K[i]] （n×1）的映射，并返回参数coef，方式如下：

用LinearRegression拟合 X[:,K[:i]] 到X[:,K[i]]的映射，并返回特征的权重系数weight_lr

用 LassoLarsIC拟合 X[:,K[:i]]×weight_lr 到X[:,K[i]]的函数映射，并返回参数weight_lasso

最终coef=weight_lr×weight_lasso

2.2更新因果关系到矩阵B

B[K[i],K[:i]]=coef

二、FGES算法

Algorithm 1: FGES算法

Input：m*n矩阵，m为样本数，n为变量数

Output：有向无环图G，G中有向边表示变量间因果关系

创建一张有向无环图G，添加n个节点，每个节点代表一个变量，每条有向边代表从原因指向结果的因果关系。

从空图两两计算添加新边后BIC得分差值，得到候选的有向边列表。<---- Algorithm 2

前向贪婪搜索，给图G添加得分差值最高的边。<---- Algorithm 4

为后向贪婪搜索准备候选有向边列表。 <---- Algorithm 9

后向贪婪搜索，给图G删除得分差值最高的边。

Algorithm 2: 从空图初始化，得到候选的有向边列表

Input：图G，G中节点集V={v_1, v_2, ..., v_n}，边集E={e_1, e_2, ..., e_k}

Output：候选的Arrow列表，元素为4元组，包含候选的有向边(x，y)及相应的候选条件集S、马尔科夫毯集T、BIC得分

记Arrow列表，元素为4元组，包含候选的有向边(x，y)及相应的候选条件集S、马尔科夫毯集T、BIC得分

For v_i, v_j in V:

计算在空图上添加v_i--v_j边后的BIC得分差值 <- --- Algorithm3

if BIC得分差值> 0:

在Arrow中分别添加从v_i到v_j，从v_j到v_i的边及BIC得分，对应候选条件集、马尔科夫毯集初始化为空集

按照BIC得分差值逆序排列Arrow

return Arrow

Algorithm 3: 计算BIC得分差值

Input：准备添加的边x--y；y候选的父集F；惩罚常数c

Output：BIC得分差值

计算y在F上回归的残差s_0，y在F∪{x}上回归的残差s_1。

计算BIC得分，BIC（F，y）= -n*ln(s^2)-c*k*ln(n)，k=|F|+1，n是样本大小，s^2是y在F上回归残差的方差，c是用于增加或减少分数惩罚的常数。

return BIC(F∪{x}, y) - BIC(F, y)

Algorithm 4: 前向贪婪搜索

input：可能的Arrow列表，元素为4元组，包含可能的有向边(x，y)及相应的候选条件集S、马尔科夫毯集T、BIC得分；总分total_score

output：加边后的图G

当Arrow不为空：

获取得分最高的关系对（x, y）,从Arrow中删除该条记录，如果x与y邻接，continue

计算候选条件集S。从图G中计算与x，y均有边，且与y为无向边的候选条件集S。检查最大度的限制max_degree，如果S成员数量大于max_degree，continue。检查图G中计算的候选条件集与Arrow中缓存的候选条件集的一致性，如果不一致，continue。

计算马尔科夫毯集T。从图G中搜索用于检查x-y条件独立性的马尔科夫毯集T <---- Algorithm 5

如果Arrow中缓存的马尔科夫毯集T不是图G中搜索的马尔科夫毯集T的子集，continue。如果可能的候选条件集S与马尔科夫毯集T的并集在图中不是全连接，continue

检查x与y之间，以S与T的并集为条件集，是否存在连通的半有向路径，如果存在，x与y不满足马尔科夫条件独立性，continue <--- Algorithm 6

在图G添加x-->y

对于马尔科夫毯集T中节点t_i，只保留t_i到y的有向边，删除y到t_i的有向边(假如存在)

总分total_score += 添加x-->y的BIC得分差值

应用meek rules对添加x-->y边后的图G重新确定方向，记录检查过的节点集P < ---- Algorithm 7

基于图G更新Arrow <---- Algorithm 8

Algorithm 5: 搜索马尔科夫毯

input：图G，外部变量x，目标变量y

output：马尔科夫毯集T

初始化马尔科夫毯集T为空集

for z in y的所有邻居节点：

if z与y有无向边：

if z与x无边：

将z加入T中

return 马尔科夫毯集T

Algorithm 6: 检查马尔科夫条件独立性

input：图G，起点x，终点y，候选条件集S

output：bool值，表示是否存在非堵塞的半有向路径

for path in x到y的所有半有向路径：

if 路径path中的节点在候选条件集S中，continue

return True

return False

Algorithm 7: 4条meek rules

input：图G，节点x

output：定向后的图G

在图G中搜索x的邻接节点集A，在A 与x构成的子图中，检查以下四种结构是否存在，按照规则确定边的方向。

R1: 如果a-->b, b---c, 且a和c不邻接，则b-->c

R2: 如果a-->b,b-->c, 且a和c邻接，则a-->c

R3: 如果a与b、c、d都邻接，b与c不邻接，且b-->d<--c，则a-->d

R4: 如果 a与b、c、d都邻接，b与c不邻接，且b-->d-->c，则a-->c

Algorithm 8: 基于图G更新Arrow（前向）

input：图G，检查过的节点集P，可能的Arrow列表

output：更新后的Arrow

for t in P：

for z in 与t的BIC得分大于0的节点：

如果z与t在图G中不邻接，在Arrow中删除z-->t的记录

计算z--t的候选候选条件集S，马尔科夫毯集T

for sub(T) in T的全组合：

将sub(T)作为t的马尔科夫毯，S为候选条件集，计算图G添加z--t后的BIC得分差值，如果差值大于0，将以sub(T)为马尔科夫毯、S为候选条件集的z--t边记录到候选Arrow列表。

Algorithm 9: 初始化后向贪婪搜索的候选Arrow列表

input：图G

output：候选Arrow列表

for x-->y in 图G:

计算x-->y的候选条件集S

for sub(S) in S的全组合：

将sub(S)作为y的马尔科夫毯，计算图G删除x-->y后BIC得分的差值 <---- Algorithm 3

如果差值大于0，将以S - sub(S)为马尔科夫毯、S为候选条件集的x-->y边记录到候选Arrow列表。

Algorithm 10: 后向贪婪搜索

input：可能的Arrow列表，元素为4元组，包含可能的有向边(x，y)及相应的候选条件集S、马尔科夫毯集T、BIC得分；总分total_score

output：删边后的图G

当Arrow不为空：

获取得分最高的关系对（x, y）,从Arrow中删除该条记录，如果x与y不邻接或存在y-->x边，continue。如果Arrow中记录的条件集S与马尔科夫毯集T的差集在图G不是完全图，continue。如果图G中计算的候选条件集与Arrow中缓存的候选条件集不一致，continue。

在图G中删除x--y，将{x,y}--T定向为{x, y}-->T

应用meek rules对删除x-->y边后的图G重新确定方向，记录检查过的节点集P < ---- Algorithm 7

总分total_score += 删除x-->y的BIC得分差值

基于图G更新Arrow <---- Algorithm 11

Algorithm 11: 基于图G更新Arrow（后向）

input：图G，检查过的节点集P，可能的Arrow列表

output：更新后的Arrow

for t in P：

for z in t的邻接集合：

if 图G中存在z-->t：

在Arrow中删除z-->t、t-->z的记录。

基于z--t扩展Arrow集 <---- Algorithm 12

elif 图G中存在z--t：

在Arrow中删除z-->t、t-->z的记录。

基于z--t扩展Arrow集 <---- Algorithm 12

基于t--z扩展Arrow集 <---- Algorithm 12

Algorithm 12：基于x-->y扩展Arrow集（后向）

input：图G，x-->y边

output：更新后的Arrow

在图G计算x-->y的候选条件集S

for sub(S) in S的全组合：

将sub(S)作为y的马尔科夫毯，计算图G删除x-->y后BIC得分的差值 <---- Algorithm

如果差值大于0，将以S - sub(S)为马尔科夫毯、S为候选条件集的x-->y边记录到候选Arrow列表。