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

2024-03-28 09:46:23 发布

您现在位置:Python中文网/ 问答频道 /正文

我正在用python实现PC算法。这种算法构造了n元高斯分布的图形模型。此图形模型基本上是有向无环图的骨架,这意味着如果结构类似:

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

在图中,那么x1与给定的x2的x3无关。更一般地说,如果A是图的邻接矩阵,并且A(i,j)=A(j,i)=0(i和j之间有一条缺失的边),那么i和j是条件独立的,由i到j的任何路径中出现的所有变量决定。出于统计和机器学习的目的,可以“学习”底层的图形模型。 如果我们对联合高斯n变量随机变量有足够的观察,我们可以使用PC算法,其工作原理如下:

^{pr2}$

该算法还计算图的分离集,该分离集由另一个算法使用,该算法从骨架和pc算法返回的分离集开始构造dag。这就是我目前所做的:

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)

该函数实现了一个统计检验,该检验使用估计偏相关的fisher z变换。我以两种方式使用此算法:

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

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


Tags: ofthetest算法lenreturnifpartial
1条回答
网友
1楼 · 发布于 2024-03-28 09:46:23

我说到点子上了。上述代码中最关键的问题是以下错误:

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是精度矩阵适当子集上的Fisher变换,它估计给定K的i和j的偏相关。正确的检验如下:

^{pr2}$

我希望有人能发现这个有用

相关问题 更多 >