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
|
function [pdag, G] = learn_struct_pdag_pc(cond_indep, n, k, varargin)
% LEARN_STRUCT_PDAG_PC Learn a partially oriented DAG (pattern) using the PC algorithm
% P = learn_struct_pdag_pc(cond_indep, n, k, ...)
%
% n is the number of nodes.
% k is an optional upper bound on the fan-in (default: n)
% cond_indep is a boolean function that will be called as follows:
% feval(cond_indep, x, y, S, ...)
% where x and y are nodes, and S is a set of nodes (positive integers),
% and ... are any optional parameters passed to this function.
%
% The output P is an adjacency matrix, in which
% P(i,j) = -1 if there is an i->j edge.
% P(i,j) = P(j,i) = 1 if there is an undirected edge i <-> j
%
% The PC algorithm does structure learning assuming all variables are observed.
% See Spirtes, Glymour and Scheines, "Causation, Prediction and Search", 1993, p117.
% This algorithm may take O(n^k) time if there are n variables and k is the max fan-in,
% but this is quicker than the Verma-Pearl IC algorithm, which is always O(n^n).
sep = cell(n,n);
ord = 0;
done = 0;
G = ones(n,n);
G=setdiag(G,0);
while ~done
done = 1;
[X,Y] = find(G);
for i=1:length(X)
x = X(i); y = Y(i);
%nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]);
nbrs = mysetdiff(neighbors(G, y), x); % bug fix by Raanan Yehezkel <raanany@ee.bgu.ac.il> 6/27/04
nbrs = unique(nbrs); % bug fix Immer Ebert ebert@tree.com 3/278/11
if length(nbrs) >= ord & G(x,y) ~= 0
done = 0;
%SS = subsets(nbrs, ord, ord); % all subsets of size ord
SS = subsets1(nbrs, ord);
for si=1:length(SS)
S = SS{si};
if feval(cond_indep, x, y, S, varargin{:})
%if isempty(S)
% fprintf('%d indep of %d ', x, y);
%else
% fprintf('%d indep of %d given ', x, y); fprintf('%d ', S);
%end
%fprintf('\n');
% diagnostic
%[CI, r] = cond_indep_fisher_z(x, y, S, varargin{:});
%fprintf(': r = %6.4f\n', r);
G(x,y) = 0;
G(y,x) = 0;
sep{x,y} = myunion(sep{x,y}, S);
sep{y,x} = myunion(sep{y,x}, S);
break; % no need to check any more subsets
end
end
end
end
ord = ord + 1;
end
% Create the minimal pattern,
% i.e., the only directed edges are V structures.
pdag = G;
[X, Y] = find(G);
% We want to generate all unique triples x,y,z
% This code generates x,y,z and z,y,x.
for i=1:length(X)
x = X(i);
y = Y(i);
Z = find(G(y,:));
Z = mysetdiff(Z, x);
for z=Z(:)'
if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x})
%fprintf('%d -> %d <- %d\n', x, y, z);
pdag(x,y) = -1; pdag(y,x) = 0;
pdag(z,y) = -1; pdag(y,z) = 0;
end
end
end
% Convert the minimal pattern to a complete one,
% i.e., every directed edge in P is compelled
% (must be directed in all Markov equivalent models),
% and every undirected edge in P is reversible.
% We use the rules of Pearl (2000) p51 (derived in Meek (1995))
old_pdag = zeros(n);
iter = 0;
while ~isequal(pdag, old_pdag)
iter = iter + 1;
old_pdag = pdag;
% rule 1
[A,B] = find(pdag==-1); % a -> b
for i=1:length(A)
a = A(i); b = B(i);
C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a
if ~isempty(C)
pdag(b,C) = -1; pdag(C,b) = 0;
%fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C);
end
end
% rule 2
[A,B] = find(pdag==1); % unoriented a-b edge
for i=1:length(A)
a = A(i); b = B(i);
if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' );
pdag(a,b) = -1; pdag(b,a) = 0;
%fprintf('rule 2: %d -> %d\n', a, b);
end
end
% rule 3
[A,B] = find(pdag==1); % a-b
for i=1:length(A)
a = A(i); b = B(i);
C = find( (pdag(a,:)==1) & (pdag(:,b)==-1)' );
% C contains nodes c s.t. a-c->ba
G2 = setdiag(G(C, C), 1);
if any(G2(:)==0) % there are 2 different non adjacent elements of C
pdag(a,b) = -1; pdag(b,a) = 0;
%fprintf('rule 3: %d -> %d\n', a, b);
end
end
end
|
