Я реализую алгоритм ПК в Python.Такой алгоритм строит графическую модель n-вариативного гауссовского распределения.Эта графическая модель в основном является каркасом ориентированного ациклического графа, что означает, что если в графе есть структура типа:
(x1)---(x2)---(x3)
, то x1 не зависит от x3, заданного x2.В более общем случае, если A является матрицей смежности графа и A (i, j) = A (j, i) = 0 (между i и j отсутствует ребро), то i и j условно независимы по всем переменнымкоторые появляются на любом пути от i до j.В статистических целях и в целях машинного обучения можно «изучить» базовую графическую модель.Если у нас достаточно наблюдений совместно гауссовой n-переменной случайной величины, мы могли бы использовать алгоритм ПК, который работает следующим образом:
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 - это альфа-параметр настройки, с помощью которого я буду проверять условную независимость, а 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 - лямбда, которая получает матрицу (обратную к ковариационной матрице), целое число 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, лежащую в основе определенного набора данных, либо потому, что знаю генеративную модель, но она не совпадает с той, которую изучает мой алгоритм.Я прекрасно знаю, что это нетривиальная задача, и я мог неправильно понять теоретическую концепцию, а также допущенную ошибку даже в тех частях кода, которые я здесь пропустил;но сначала я хотел бы знать (от кого-то, кто более опытен, чем я), если тест, который я написал, является правильным, а также, если есть библиотечные функции, которые выполняют такие тесты, я попытался выполнить поиск, но не смог найтилюбая подходящая функция.