D-Separation和PC算法

简介

综合转载于:

D-Separation

很多的机器学习模型都可以用概率的角度去解释(可以看MLAPP和PRML这两本书),其中一类重要的模型就是概率图模型,而是概率图模型的灵魂就是模型变量间的条件独立性。 因为有了独立性,才有了各种不同的概率图模型,比如LDA,HMM等等模型。那么概率图中,变量间的独立性是怎么体现的呢?D-分离准则就是一种简单的技巧去判断一个概率图中的独立性的。

简单的说,如果在概率图模型中,比如说X,Y两个节点没有边,使得X和Y之间肯定存在某种独立性,这种独立性可以是在某个子集Z的条件下使得他们独立,也可能是他们两个本身就是独立的,这时我们称X,Y之间是D-分离的。现在先给出D-分离的准则:

定义: 当路径p被结点集Z,d-分离(或被blocked掉)时,当且仅当以下条件成立: 1. 若p包含形式如下的链i->m->j 或i<-m->j,则结点m在集合Z中。 2. 若p中包含collider(碰撞点)i->m<-j,则结点m不在Z中且m的子代也不在Z中。

更进一步说,如果Z将X和Y d-separate,当且仅当Z将X,Y之间的每一条路径都block掉。

接下来逐步的介绍上述的准则,我们可以拆分成3个规则来考虑:

首先,这里我们先说明一下path,我们说两个结点之间的path的时候是不管他们之间边的方向的。

没有条件集的独立性

规则1如果x到y的任一path(路径)都经过collider(碰撞点),则x和y独立。注意,这里的路径是忽略边的方向的,而碰撞点是指有多个箭头指向的它的节点,即类似于下图的s->t<-u。

现在我们看看这个图变量之间的独立性是怎样的。先考虑x和t的路径:x-r-s-t。这条路径中并不存在碰撞点,所以x和t不独立,同样地,t和y、u、v都不独立。然而,对于变量x到y而言,x,y之间的路径必然会经过碰撞点t,所以x和y是独立,同样地,在碰撞点两侧的变量都是独立的,比如x和v,s和u,r和u也是独立的。

一般的条件独立

那么如果一条路径中没有碰撞点,怎么才能让他们独立呢?我们还能使用条件独立的性质,只要条件集Z能够将一条路径block掉,那么就可以让两个变量独立。注意的是,当碰撞点或碰撞点的子代出现在条件集的时候要小心,很有可能会导致不同的结果,这个问题我们留到规则3。现在我们先假设条件集中不存在碰撞点。

规则2当x到y的之间的任一路径都经过Z中的节点,且Z并不包含碰撞点或碰撞点的子代,则x和y独立。

如图,设结点集Z={r, v}(图中画圈的节点),根据规则2,x和s是条件独立的,因为x和s之间的路径被block掉了,同样地,u和y也是条件独立的。

当collider作为条件集

现在我们先看看collider成为条件集的时候会发生什么。我们考虑一个例子

这里显然IQ和难度是独立的,他们共同决定成绩的高低,于是成绩是一个collider,而不同成绩会导致不同排名,所以排名是成绩的子代。现在如果我们已经知道一个人的成绩是多少,那么神奇的事情发生了,如果我们知道一个人IQ很低,但成绩很高,那么我们可以推断出难度肯定很低;同理我们知道成绩很高,但难度很大,那么可以推断出IQ很高。原本IQ和难度两个独立的变量在成绩的条件下变得不独立的。

同理,因为ABCD等级是由不同成绩划分出来的,也包含了成绩的信息,所以我们知道一个人的等级同样可以推断出一个人的IQ和难度。因此当条件集是collider的时候,会反而使得两个独立的变量变得不独立。这就是我们判断条件独立性的时候要小心的地方。接下来是规则3,

规则3当碰撞点或碰撞点的子孙节点为集合Z的成员时,该碰撞点不再截断路径。

设结点集Z={r, p}(图中画圈的节点),根据规则3,在给定Z的情况下,s和y不独立,因为t的孩子节点p在集合Z中,所以t没有办法像规则1一样截断路径s-t-u-v-y,与此相反,条件集p使得s和u变得不独立了。然而在这里,x和u是独立的,虽然t不能截断它,但是r可以截断它(根据规则2)。

PC 算法

概率图模型

对于一个 $K$ 维随机向量 $X=\left[X_1, X_2, \ldots, X_K\right]^{\top}$ ,一般难以直接建模。因为如果每个变量为离散变量并有 $M$ 个可能取值,在不作任何独立性假设的前提下,需要 $M^K-1$ 个参数才能表示其概率分布,参数数量会非常庞 大。

一种减少参数数量的方法是独立性假设。把 $X$ 的联合概率分解为 $K$ 个条件概率的乘积:

$x$ 为随机向量 $X$ 的取值。可以看到,如果某些变量之间存在条件独立,参数数量就可以大幅咸少。
因此,概率图模型(Probabilistic Graphical Model,PGM) 用图结构来描述多元随机变量之间的条件独立关系, 从而为研究高维空间中的概率模型带来了很大的便婕。
概率图模型中,每个节点表示一个 (或一组) 随机变量,边表示文些随机变量之间的概率依赖关系。常见的概率图 模型可以分为有向图模型和无向图模型。

  • 有向图模型 (Directed Graphical Model),也称为贝叶斯网络 (Bayesian Network) 或信念网络 (Belief Network) ,使用有向无不图 (Directed Acyclic Graph,DAG) 来描述变量之间的关系。如果两个节点之间 有连边,表示文两个变量之间有因果关系,即不存在其他变量使得这两个变量条件独立。
  • 无向图模型,也称为马尔可夫陇机场(Markov Random Field,MRF),使用无向图来描述变量之间的关 系。两个节点之间有连边代表这两个变量之间有概率依赖关系,但不一定是因果关系。

本文只讨论有向图模型,即贝叶斯网络。

贝叶斯网络

定义

有向无环图 $G$ 中,每个节点对应 $K$ 维随机向量 $X$ 中的一个变量,有向边 $e_{i j}$ 表示随机变量 $X_i$ 和 $X_j$ 之间具有 因果关系,父节点 $X_i$ 是『因』,子节点 $X_j$ 是『果』,显然这两个点之间一定非条件独立。
令 $X_{\pi_k}$ 为变量 $X_k$ 的所有父节点变量集合, $P\left(X_k \mid X_{\pi_k}\right)$ 表示每个随机变量的局部绦件㪜率分布 (Local Conditional Probability Distribution) 。
如果 $X$ 的联合概率分布可以分解为每个随机变量 $X_k$ 的局部条件概率的连乘形式,即:

那么称 $(G, X)$ 构成了一个贝叶斯网络。

局部马尔可夫性质

每个随机变量在给定父节点的情况下,条件独立于它的非后代节点:

其中 $Z$ 为 $X_k$ 的非后代节点。

基本问题

  • 学习问题

    • 结构学习:那么怎样才可以获得这个神奇的有向无环图呢,这就是结构学习问题。即学习出最优网络结构,也就是各节点之间的依赖关系。主流的结构学习方法主要可以分为:

      • 基于评分搜索的方法:利用搜索算法和评分函数,对每一个搜索到的网络结构进行评分,最终搜索出评分最高的网络结构。搜索算法的复杂度和精确度直接决定了学习算法的搜索效率,评分函数的优劣也直接决定了算法的计算复杂度和精确度。所以选择合理的优化搜索算法和评分函数是该类方法的核心问题。

        该类方法容易陷入局部最优解而无法达到全局最优,并且结构空间的⼤⼩随节点的增加呈指数增加(空间复杂度)。

      • 基于依赖统计分析的方法:分析变量间的依赖关系,在依赖性较大的两节点之间添加连接边,得到无向图。然后根据包含关系等方式确定无向图中边的方向,得到最终的有向无环图。本文要讨论的 PC 算法就是这类方法中(比较古老)的一种。

        该类方法能获得全局最优解,但随着节点的增加,算法的时间复杂度会增加得很快;并且它不能区分同属于一个马尔可夫等价类的图,这一点后面会讲到。

    • 参数学习:在给定网络结构时,确定网络参数,即参数估计问题:

      • 不含隐变量:如果图模型中不含隐变量(latent variable),即所有变量都是可观测的,那么网络参数一般可以直接通过最大似然来进行估计。

      • 含隐变量:有些时候 X 中的变量有很复杂的依赖关系,这时通常会引入隐变量 z 来简化模型。如果图模型中包含隐变量,即有部分变量是不可观测的,这时就需要用 EM 算法(Expectation Maximum,期望最大化算法)来进行参数估计。

        如果 EM 算法中的后验是 intractable 的,那么又需要用变分推断(Variational Inference)来寻找一个简单分布来近似后验。

        而在深度学习大行其道的今天,你可能会想到用神经网络去拟合这个后验不就完事儿了,是的这就是变分自编码器(Variational Auto-Encoder,VAE)的思想,去学它吧朋友。

  • 推断问题:在已知部分变量时,计算其他变量的条件概率分布

PC 算法

好的现在讲主题了,用 PC 算法来学习出贝叶斯网络的结构。如上文所述,PC 算法会先确定节点间的依赖关系(但不确定方向),即先生成一个无向图,然后再确定依赖方向,把无向图扩展为完全部分有向无环图(Completed Partially Directed Acyclic Graph,CPDAG)。

依赖关系确立

设 $V$ 是输入点集,有以下步㵵:

  • 在 $V$ 上生成完全无向图 $G$
  • 对于 $G$ 中的两个相邻点 $i, j$ ,如果 $i$ 和 $j$ 能在给定节点 $k$ 时条件独立,则删除 $i$ 和 $j$ 之间的边
    这样会得到一个无向图,图中的无向边表示它连接的两个节点之间有依赖 (因果) 关系,这样的无向图叫骨喿 (skeleton)。PC 算法把上述过程转化为了 $d$ 分隔 (d-separation) 问题。

d 分隔

d 分隔的定义

节点集合 Od 分隔节点 i 与节点 j,当且仅当:

给定 O 时,ij 之间不存在有效路径(active path),即 ijO 下条件独立(记作 ijO)。

用 $O(i, j)$ 表示能够 $d$ 分隔 $i$ 和 $j$ 的点集,用 $a d j(G, x)$ 表示图 $G$ 中节点 $x$ 的相邻点集,那么 PC 算法检验条件独立性的具体流程为:

简单总结一下:

  • $\ell=1$
  • repeat
    • for 每个相邻点对 $(i, j)$
      • for $\operatorname{adj}(G, i) \backslash\{j\}$ 或 $\operatorname{adj}(G, i) \backslash\{i\}$ 的所有可能的节点数为 $\ell$ 的子集 $K$
        • 测试 $K$ 能否 $d$ 分隔 $(i, j)$
        • 如果能,说明 $i$ 和 $j$ 之间不存在有效的依赖关系,所以删除边 $i-j$ ,并将这个点集加入 $O(i, j)$ 和 $O(j, i)$, break
    • $\ell=\ell+1$
  • until 当前图中的所有的邻接点集都小于 $\ell$

Fisher Z Test

为了判断 $d$ 分隔,我们需要对任意两个节点进行条件独立性检验, PC 算法采用了 Fisher Z Test 作为条件独立性 检验方法。实际上 Fisher Z Test 是一种相关性检验方法,但 PC 算法认为这一堆随机变量整体上服从多元高斯分 布,这时变量条件独立与变量之间的偏相关系数为 0 等价(多元高斯分布的基本特性,证明过程可以参考 Steffen L. Lauritzen 的课件第 4.2.1 节),所以可以用 Fisher Z Test 进行条件独立性检验。

偏相关系数指校正其它变量后某一变量与另一变量的相关关系,校正的意思可以理解为假定其它变量都取值为均 数。任意两个变量 $i, j$ 的 $h$ 阶(排除其他 $h$ 个变量的影响后, $h<=k-2$ ) 偏相关系数为:

为了判断 $\rho$ 是否为 0 ,需要将 $\rho$ 通过 Fisher Z 变换 $[6]$ 转换成正态分布:

定义零假设和对立假设:

  • 零假设: $H_0(i, j \mid K): \rho_{i, j \mid K} \neq 0$
  • 对立假设: $H_1(i, j \mid K): \rho_{i, j \mid K}=0$

然后给定一个显著性水平 $\alpha \in(0,1)$ ,那么 (双侧) 检验的规则为,如果有:

其中 $\Phi(\cdot)$ 为 $\mathcal{N}(0,1)$ 的累积分布函数,则拒绝零假设, $i, k$ 关于 $K$ 条件独立。所以将上面伪代码的第 11 行替 换为 if $\sqrt{n-|K|-3} \mid Z(i, j \mid K) \leq \Phi^{-1}(1-\alpha / 2)$ 。

依赖关系方向确立

经过上一个阶段,我们得到了一个无向图。现在我们要利用 $d$ 分隔的原理来确定图中边的依赖方向,把骨架扩展 为 $\mathrm{DAG}$ 。

对于任意三个以有效依赖关系边相连的节点 $X-Z-Y$ ,其依赖关系必为下图的四种关系之一:

$d$ 分隔的结论为: 对于有向无环图 $E$ ,有两个节点 $X, Y$ 和一个点集 $O$ ,为了判断 $X$ 和 $Y$ 是否关于 $O$ 条件独 立,考虑 $E$ 中所有 $X$ 和 $Y$ 之间的无向路径,对于其中一条路径,如果它满足以下两个条件中的任意一条,则称 这条路径是阻塞的:

  • 路径中存在某个节点 $Z$ 是 head-to-tial (图中情兄 $\mathrm{a}$, b) 或 tail-to-tail 节点 (图中情况 c) , 且 $Z$ 包含在 $O$ 中
  • 路径中存在某个节点 $Z$ 是 head-to-head 节点 (图中情兄 d),且 $Z$ 没有被包含在 $O$ 中

如果 $X, Y$ 间所有的路径都是阻塞的,那么 $X, Y$ 关于 $O$ 条件独立;否则, $X, Y$ 不关于 $O$ 条件独立。

而我们已经记录了 $d$ 分隔 $X$ 和 $Y$ 的点集 $O$ ,因此我们可以由 $d$ 分隔的结论反推出贝叶斯网络中边的方向,方向 的判断方法可以转换成以下三条规则:

  • 规则 1: 如果存在 $X \rightarrow Y-Z$ ,把 $Y-Z$ 变为 $Y \rightarrow Z$
  • 规则 2: 如果存在 $X \rightarrow Z \rightarrow Y$ ,把 $X-Y$ 变为 $X \rightarrow Y$
  • 规则 3: 如果存在 $X-Z_1 \rightarrow Y, X-Z_2 \rightarrow Y$, 且 $Z_1, Z_2$ 不相邻,把 $X-Y$ 变为 $X \rightarrow Y$

实际上还可以推出一个规则 4 :

  • 规则 4: 如果存在 $X-Z_1 \rightarrow Z_2$ 和 $Z_1 \rightarrow Z_2 \rightarrow Y$ ,且 $Z_1, Z_2$ 不相邻,把 $X-Y$ 变为 $X \rightarrow Y$

但很显然这种情兄是矛盾的,不可能存在,所以不用考虑。

总结一下:

这样我们就可以得到一个完全部分有向无环图。

马尔科夫等价类

这部分我就没看懂了。

很明显,完全部分有向无环图(CPDAG)跟有向无环图看上去就不一样。首先来看什么是部分有向无环图(Partially Directed Acyclic Graph,PDAG):

部分有向无环图

假设 $G=(V, E)$ 是一个图,若氻集 $E$ 中包含有向边和无向边,且不存在有向环,则称 $G$ 是一个部分有 向无环图。

完全部分有向无环图指:

完全部分有向无环图

假设 $G=(V, E)$ 是一个部分有向无环图,若 $E$ 中的有向边都是不可逆的,并且 $E$ 中的无向边都是可 逆的,则称 $G$ 是一个完全部分有向无环图。

关于可逆和不可逆:

可逆 / 不可逆

对于有向无环图 $G=(V, E)$ 中的任意有向边 $V_i \rightarrow V_j \in E$ ,如果存在图 $G^{\prime}=\left(V, E^{\prime}\right)$ 与 $G$ 等价, 且 $V_j \rightarrow V_i \in E^{\prime}$ ,则称有向边 $V_i \rightarrow V_j$ 在 $G$ 中是可逆的,否则是不可逆的。
同理,对任意无向边 $V_i-V_j \in E$ ,若存在 $G_1=\left(V, E_1\right) 、 G_2=\left(V, E_2\right)$ 均与 $G$ 等价,且 $V_i \rightarrow$ $V_j \in E_1 、 V_j \rightarrow V_i \in E_2$ ,则称无向边 $V_i-V_j$ 在 $G$ 中是可逆的,否则是不可逆的。

换句话说用 PC 算法得到的图是含有无向边的。这是因为依据 $d$ 分隔确定的条件独立性所构造的网络结构不具有 唯一性,它们只是直实的贝叶斯网络的马尔科夫等价类(Markov Equivalence Class):

马尔科夫等价类

有向无环图 $G_1=\left(V, E_1\right)$ 和 $G_2=\left(V, E_2\right)$ 有相同的顶点集合和骨架, $V$ 为顶点集合, $E_1$ 和 $E_2$ 为 边的集合。
对于任意的不相交的顶点集合 $A, B, C \in V$ ,如果满足 $A, B$ 在 $G_1$ 和 $G_2$ 中都被 $C$ 所 $d$ 分隔(也叫有 相同的 $V$ 结构) ,则称图 $G_1$ 和 $G_2$ 是马尔科夫等价的。

举个栗子:

上图 $G_1$ 和 $G_2$ 是马尔科夫等价类,它们左上角的那条有向边方向并不相同,这时 PC 算法就无法判断这条边的方 向了,只能输出无向边,即 $G_3$ 。

所以,严格来说,PC 算法以及大多数基于依赖统计分析的贝叶斯网络结构学习算法,得到的都只是一个 CPDAG (依然有无向边),而不是真正意义上的贝叶斯网络 (有向无环图)。

参考

  1. An Algorithm for Fast Recovery of Sparse Causal Graphs. Peter Spirtes and Clark Glymour. Social Science Computer Review 1991.
  2. d-Separation: From Theorems to Algorithms. Dan Geiger, et al. UAI 1989.
  3. Estimating High-Dimensional Directed Acyclic Graphs with the PC-Algorithm. Markus Kalisch and Peter Buhlmann. JMLR 2007.
  4. Frequency Distribution of the Values of the Correlation Coefficient in Samples from an Indefinitely Large Population. R. A. Fisher. Biometrika 1915.
  5. Elements of Graphical Models. Steffen L. Lauritzen. 2011.
  6. Wikipedia: Fisher transformation

代码

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
from itertools import combinations
from scipy.stats import norm
import numpy as np
import math
from typing import List

def get_neighbors(G, x: int, y: int):
return [i for i in range(len(G)) if G[x][i] == True and i != y]

def gauss_ci_test(suff_stat, x: int, y: int, K: List[int], cut_at: float = 0.9999999):
"""条件独立性检验"""
C = suff_stat["C"]
n = suff_stat["n"]

# ------ 偏相关系数 ------
if len(K) == 0: # K 为空
r = C[x, y]

elif len(K) == 1: # K 中只有一个点,即一阶偏相关系数
k = K[0]
r = (C[x, y] - C[x, k] * C[y, k]) / math.sqrt((1 - C[y, k] ** 2) * (1 - C[x, k] ** 2))

else: # 其实我没太明白这里是怎么求的,但 R 语言的 pcalg 包就是这样写的
m = C[np.ix_([x] + [y] + K, [x] + [y] + K)]
p = np.linalg.pinv(m)
r = -p[0, 1] / math.sqrt(abs(p[0, 0] * p[1, 1]))

r = min(cut_at, max(-cut_at, r))

# Fisher's z-transform
z = 0.5 * math.log1p((2 * r) / (1 - r))
z_standard = z * math.sqrt(n - len(K) - 3)

# Φ^{-1}(1-α/2)
p_value = 2 * (1 - norm.cdf(abs(z_standard)))

return p_value

def skeleton(suff_stat, alpha: float):
n_nodes = suff_stat["C"].shape[0]

# 分离集
O = [[[] for _ in range(n_nodes)] for _ in range(n_nodes)]

# 完全无向图
G = [[i != j for i in range(n_nodes)] for j in range(n_nodes)]

# 点对(不包括 i -- i)
pairs = [(i, (n_nodes - j - 1)) for i in range(n_nodes) for j in range(n_nodes - i - 1)]

done = False
l = 0 # 节点数为 l 的子集

while done != True and any(G):
done = True

# 遍历每个相邻点对
for x, y in pairs:
if G[x][y] == True:
neighbors = get_neighbors(G, x, y) # adj(C,x) \ {y}

if len(neighbors) >= l: # |adj(C, x) \ {y}| > l
done = False

# |adj(C, x) \ {y}| = l
for K in set(combinations(neighbors, l)):
# 节点 x, y 是否被节点数为 l 的子集 K d-seperation
# 条件独立性检验,返回 p-value
p_value = gauss_ci_test(suff_stat, x, y, list(K))

# 条件独立
if p_value >= alpha:
G[x][y] = G[y][x] = False # 去掉边 x -- y
O[x][y] = O[y][x] = list(K) # 把 K 加入分离集 O
break

l += 1

return np.asarray(G, dtype=int), O

def extend_cpdag(G, O):
n_nodes = G.shape[0]

def rule1(g):
"""Rule 1: 如果存在链 i -> j - k ,且 i, k 不相邻,则变为 i -> j -> k"""
pairs = [(i, j) for i in range(n_nodes) for j in range(n_nodes) if g[i][j] == 1 and g[j][i] == 0] # 所有 i - j 点对

for i, j in pairs:
all_k = [k for k in range(n_nodes) if (g[j][k] == 1 and g[k][j] == 1) and (g[i][k] == 0 and g[k][i] == 0)]

if len(all_k) > 0:
g[j][all_k] = 1
g[all_k][j] = 0

return g

def rule2(g):
"""Rule 2: 如果存在链 i -> k -> j ,则把 i - j 变为 i -> j"""
pairs = [(i, j) for i in range(n_nodes) for j in range(n_nodes) if g[i][j] == 1 and g[j][i] == 1] # 所有 i - j 点对

for i, j in pairs:
all_k = [k for k in range(n_nodes) if (g[i][k] == 1 and g[k][i] == 0) and (g[k][j] == 1 and g[j][k] == 0)]

if len(all_k) > 0:
g[i][j] = 1
g[j][1] = 0

return g

def rule3(g):
"""Rule 3: 如果存在 i - k1 -> j 和 i - k2 -> j ,且 k1, k2 不相邻,则把 i - j 变为 i -> j"""
pairs = [(i, j) for i in range(n_nodes) for j in range(n_nodes) if g[i][j] == 1 and g[j][i] == 1] # 所有 i - j 点对

for i, j in pairs:
all_k = [k for k in range(n_nodes) if (g[i][k] == 1 and g[k][i] == 1) and (g[k][j] == 1 and g[j][k] == 0)]

if len(all_k) >= 2:
for k1, k2 in combinations(all_k, 2):
if g[k1][k2] == 0 and g[k2][k1] == 0:
g[i][j] = 1
g[j][i] = 0
break

return g

# Rule 4: 如果存在链 i - k -> l 和 k -> l -> j,且 k 和 l 不相邻,把 i - j 改为 i -> j
# 显然,这种情况不可能存在,所以不需要考虑 rule4

# 相邻点对
pairs = [(i, j) for i in range(n_nodes) for j in range(n_nodes) if G[i][j] == 1]

# 把 x - y - z 变为 x -> y <- z
for x, y in sorted(pairs, key=lambda x:(x[1], x[0])):
all_z = [z for z in range(n_nodes) if G[y][z] == 1 and z != x]

for z in all_z:
if G[x][z] == 0 and y not in O[x][z]:
G[x][y] = G[z][y] = 1
G[y][x] = G[y][z] = 0

# Orientation rule 1 - rule 3
old_G = np.zeros((n_nodes, n_nodes))

while not np.array_equal(old_G, G):
old_G = G.copy()

G = rule1(G)
G = rule2(G)
G = rule3(G)

return np.array(G)

def pc(suff_stat, alpha: float = 0.05, verbose: bool = False):
G, O = skeleton(suff_stat, alpha) # 骨架
cpdag = extend_cpdag(G, O) # 扩展为 CPDAG

if verbose:
print(cpdag)

return cpdag
一分一毛,也是心意。