我正在用Python实现PC算法。这种算法构建一个n变量高斯分布的图模型。这个图模型基本上是一个有向无环图的骨架，这意味着如果图中存在像这样的结构：
(x1)---(x2)---(x3)那么给定x2时，x1和x3是独立的。更一般地说，如果A是图的邻接矩阵，并且A(i,j)=A(j,i) = 0（i和j之间缺少边），那么i和j在从i到j的任何路径上出现的所有变量条件下是条件独立的。为了统计和机器学习的目的，可以“学习”底层的图模型。如果我们有足够的联合高斯n变量随机变量的观测，我们可以使用如下工作的PC算法：
给定n为观测到的变量数量，初始化图为G=K(n)对于节点的每一对i,j： 如果存在从i到j的边e： 查找i的邻居 如果j在i的邻居中则从邻居集合中移除j 将邻居集合称为k 测试给定集合k时i和j是否独立，如果为TRUE： 从i到j移除边e该算法还计算图的分离集，这些分离集被另一个从骨架和PC算法返回的分离集开始构建DAG的算法使用。这是我迄今为止所做的：
def _core_pc_algorithm(a,sigma_inverse):l = 0N = 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 + 1return (act_g,sep_set)a是用于测试条件独立性的调节参数alpha，sigma_inverse是采样观测的协方差矩阵的逆矩阵。此外，我的测试是：
def test(sigma_inverse,z,i,j,l,a,k): def erfinv(x): #用于近似高斯累积密度函数的逆 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): #计算给定一个条件变量K时i和j的部分相关性 part_corr_coeff_ij = z(sigma_inverse,i,j) #这给出i和j的部分相关系数 part_corr_coeff_iK = z(sigma_inverse,i,K) #这给出i和k的部分相关系数 part_corr_coeff_jK = z(sigma_inverse,j,K) #这给出j和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)) #这给出给定K时i和j的部分相关系数 return part_corr_coeff_ijK == 0 #给定K时i独立于j，如果部分相关系数(i,k)|K == 0（在联合高斯假设下）[可以检查绝对值是否小于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独立于j <=> 部分相关系数(i,j) == 0（在联合高斯假设下）[可以检查绝对值是否小于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]) #假设给定Y,Z时IJ独立 <=> 给定Y和给定Z时IJ独立 else: #我必须使用z-fisher函数进行独立性测试 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，要么是因为我知道生成模型但它与我的算法学习的不符。我完全知道这是一项非平凡的任务，我可能在理论概念上误解了，也可能在我在此省略的代码部分犯了错误；但首先我想知道（从比我更有经验的人那里），我写的测试是否正确，以及是否有执行这种测试的库函数，我尝试搜索但找不到任何合适的函数。
回答：
我直入主题。上述代码中最关键的问题涉及以下错误：
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的部分相关性。正确的测试如下：
if len(K) == 0: #CM是相关矩阵，我们没有条件变量（K的长度为0） r = CM[i, j] #r是i和j的部分相关性 elif len(K) == 1: #我们有一个条件变量，与之前的版本没有太大不同，除了我无需计算相关矩阵，因为我从它开始，并且pandas在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是给定K时i和j的部分相关性else: #超过一个条件变量 CM_SUBSET = CM[np.ix_([i]+[j]+K, [i]+[j]+K)] #我正在寻找的相关矩阵的子集 PM_SUBSET = np.linalg.pinv(CM_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)) #使用Fisher变换估计部分相关性return 2 * (1 - norm.cdf(abs(res))) #获得p值希望有人能发现这有帮助
