作为 PC 算法的一部分，测试 python 中的条件独立性

Question

我正在 python 中实现 PC 算法。这种算法构建了一个 n 变量高斯分布的图形模型。这个图形模型基本上是一个有向无环图的骨架，这意味着如果一个结构像：

(x1)---(x2)---(x3)

在图中，则x1独立于给定x2的x3。更一般地，如果 A 是图的邻接矩阵并且 A(i,j)=A(j,i) = 0（i 和 j 之间缺少一条边），则 i 和 j 是条件独立的，由所有变量出现在从 i 到 j 的任何路径中。出于统计和机器学习目的，可以 "learn" 底层图形模型。如果我们对联合高斯 n 变量随机变量有足够的观察，我们可以使用如下工作的 PC 算法：

given n as the number of variables observed, initialize the graph as G=K(n) 
for each pair i,j of nodes:
    if exists an edge e from i to j:
        look for the neighbours of i
        if j is in neighbours of i then remove j from the set of neighbours
        call the set of neighbours k
        TEST if i and j are independent given the set k, if TRUE:
             remove the edge e from i to j

该算法还计算图的分离集，该分离集被另一种算法使用，该算法从骨架开始构造 dag 和 pc 算法返回的分离集。这是我到目前为止所做的：

def _core_pc_algorithm(a,sigma_inverse):
l = 0
N = len(sigma_inverse[0])
n = range(N)
sep_set = [ [set() for i in n] for j in n]
act_g = complete(N)
z = lambda m,i,j : -m[i][j]/((m[i][i]*m[j][j])**0.5)
while l<N:
    for (i,j) in itertools.permutations(n,2):
        adjacents_of_i = adj(i,act_g)
        if j not in adjacents_of_i:
            continue
        else:
            adjacents_of_i.remove(j)
        if len(adjacents_of_i) >=l:
            for k in itertools.combinations(adjacents_of_i,l):
                if N-len(k)-3 < 0:
                    return (act_g,sep_set)
                if test(sigma_inverse,z,i,j,l,a,k):
                    act_g[i][j] = 0
                    act_g[j][i] = 0
                    sep_set[i][j] |= set(k)
                    sep_set[j][i] |= set(k)
    l = l + 1
return (act_g,sep_set)

a 是调整参数 alpha，我将使用它来测试条件独立性，sigma_inverse 是采样观察的协方差矩阵的逆矩阵。而且，我的测试是：

def test(sigma_inverse,z,i,j,l,a,k):
    def erfinv(x): #used to approximate the inverse of a gaussian cumulative density function
        sgn = 1
        a = 0.147
        PI = numpy.pi
        if x<0:
            sgn = -1
        temp = 2/(PI*a) + numpy.log(1-x**2)/2
        add_1 = temp**2
        add_2 = numpy.log(1-x**2)/a
        add_3 = temp
        rt1 = (add_1-add_2)**0.5
        rtarg = rt1 - add_3
        return sgn*(rtarg**0.5)
    def indep_test_ijK(K): #compute partial correlation of i and j given ONE conditioning variable K
        part_corr_coeff_ij = z(sigma_inverse,i,j) #this gives the partial correlation coefficient of i and j
        part_corr_coeff_iK = z(sigma_inverse,i,K) #this gives the partial correlation coefficient of i and k
        part_corr_coeff_jK = z(sigma_inverse,j,K) #this gives the partial correlation coefficient of j and k
        part_corr_coeff_ijK = (part_corr_coeff_ij - part_corr_coeff_iK*part_corr_coeff_jK)/((((1-part_corr_coeff_iK**2))**0.5) * (((1-part_corr_coeff_jK**2))**0.5)) #this gives the partial correlation coefficient of i and j given K
        return part_corr_coeff_ijK == 0 #i independent from j given K if partial_correlation(i,k)|K == 0 (under jointly gaussian assumption) [could check if abs is < alpha?]
    def indep_test():
        n = len(sigma_inverse[0])    
        phi = lambda p : (2**0.5)*erfinv(2*p-1)
        root = (n-len(k)-3)**0.5
        return root*abs(z(sigma_inverse,i,j)) <= phi(1-a/2)
if l == 0:
    return z(sigma_inverse,i,j) == 0 #i independent from j <=> partial_correlation(i,j) == 0 (under jointly gaussian assumption) [could check if abs is < alpha?]
elif l == 1:
    return indep_test_ijK(k[0])
elif l == 2:
    return indep_test_ijK(k[0]) and indep_test_ijK(k[1]) #ASSUMING THAT IJ ARE INDEPENDENT GIVEN Y,Z <=> IJ INDEPENDENT GIVEN Y AND IJ INDEPENDENT GIVEN Z  
else: #i have to use the independent test with the z-fisher function
    return indep_test()

其中 z 是一个 lambda，它接收一个矩阵（协方差矩阵的逆矩阵）、一个整数 i、一个整数 j 并根据以下规则计算 i 和 j 的偏相关给定所有其余变量（我在老师的幻灯片中读到的）：

corr(i,j)|REST = -var^-1(i,j)/sqrt(var^-1(i,i)*var^-1(j,j))

这个应用程序的主要核心是 indep_test() 函数：

    def indep_test():
        n = len(sigma_inverse[0])    
        phi = lambda p : (2**0.5)*erfinv(2*p-1)
        root = (n-len(k)-3)**0.5
        return root*abs(z(sigma_inverse,i,j)) <= phi(1-a/2)

这个函数实现了一个统计检验，它使用估计偏相关的费希尔 z 变换。我以两种方式使用这个算法：

从线性回归模型生成数据并将学习到的 DAG 与预期的进行比较
读取数据集并学习底层 DAG

在这两种情况下，我并不总能得到正确的结果，要么是因为我知道某个数据集下的 DAG，要么是因为我知道生成模型但它与我的算法学习的模型不一致。我完全知道这是一项非常重要的任务，我可能误解了理论概念，甚至在我在这里省略的部分代码中也犯了错误；但首先我想知道（从比我更有经验的人那里），如果我写的测试是正确的，还有是否有执行这种测试的库函数，我尝试搜索但我找不到任何合适的函数。

Answer 1

进入正题。上面代码中最关键的问题，关于以下错误：

sqrt(n-len(k)-3)*abs(z(sigma_inverse[i][j])) <= phi(1-alpha/2)

我弄错了 n 的平均值，它不是精度矩阵的大小，而是多变量观测值的总数（在我的例子中，是 10000 而不是 5）。另一个错误的假设是 z(sigma_inverse[i][j]) 必须在给定所有其他条件的情况下提供 i 和 j 的偏相关。这是不正确的，z 是在给定 K 的情况下估计 i 和 j 的偏相关的精度矩阵的适当子集上的 Fisher 变换。正确的测试如下：

if len(K) == 0: #CM is the correlation matrix, we have no variables conditioning (K has 0 length)
    r = CM[i, j] #r is the partial correlation of i and j 
elif len(K) == 1: #we have one variable conditioning, not very different from the previous version except for the fact that i have not to compute the correlations matrix since i start from it, and pandas provide such a feature on a DataFrame
    r = (CM[i, j] - CM[i, K] * CM[j, K]) / math.sqrt((1 - math.pow(CM[j, K], 2)) * (1 - math.pow(CM[i, K], 2))) #r is the partial correlation of i and j given K
else: #more than one conditioning variable
    CM_SUBSET = CM[np.ix_([i]+[j]+K, [i]+[j]+K)] #subset of the correlation matrix i'm looking for
    PM_SUBSET = np.linalg.pinv(CM_SUBSET) #constructing the precision matrix of the given subset
    r = -1 * PM_SUBSET[0, 1] / math.sqrt(abs(PM_SUBSET[0, 0] * PM_SUBSET[1, 1]))
r = min(0.999999, max(-0.999999,r)) 
res = math.sqrt(n - len(K) - 3) * 0.5 * math.log1p((2*r)/(1-r)) #estimating partial correlation with fisher's transofrmation
return 2 * (1 - norm.cdf(abs(res))) #obtaining p-value

我希望有人能觉得这有用

作为 PC 算法的一部分，测试 python 中的条件独立性

Test for conditional independence in python as part of the PC algorithm

python

machine-learning

statistics

directed-acyclic-graphs

bayesian-networks