wolffd@0
|
1 function [pdag, G] = learn_struct_pdag_ic_star(cond_indep, n, k, varargin)
|
wolffd@0
|
2 % LEARN_STRUCT_PDAG_IC_STAR Learn a partially oriented DAG (pattern) with latent
|
wolffd@0
|
3 % variables using the IC* algorithm
|
wolffd@0
|
4 % P = learn_struct_pdag_ic_star(cond_indep, n, k, ...)
|
wolffd@0
|
5 %
|
wolffd@0
|
6 % n is the number of nodes.
|
wolffd@0
|
7 % k is an optional upper bound on the fan-in (default: n)
|
wolffd@0
|
8 % cond_indep is a boolean function that will be called as follows:
|
wolffd@0
|
9 % feval(cond_indep, x, y, S, ...)
|
wolffd@0
|
10 % where x and y are nodes, and S is a set of nodes (positive integers),
|
wolffd@0
|
11 % and ... are any optional parameters passed to this function.
|
wolffd@0
|
12 %
|
wolffd@0
|
13 % The output P is an adjacency matrix, in which
|
wolffd@0
|
14 % P(i,j) = -1 if there is either a latent variable L such that i <-L-> j
|
wolffd@0
|
15 % OR there is a directed edge from i->j.
|
wolffd@0
|
16 % P(i,j) = -2 if there is a marked directed i-*>j edge.
|
wolffd@0
|
17 % P(i,j) = P(j,i) = 1 if there is and undirected edge i--j
|
wolffd@0
|
18 % P(i,j) = P(j,i) = 2 if there is a latent variable L such that i<-L->j.
|
wolffd@0
|
19 %
|
wolffd@0
|
20 % The IC* algorithm learns a latent structure associated with a set of observed
|
wolffd@0
|
21 % variables.
|
wolffd@0
|
22 % The latent structure revealed is the projection in which every latent variable is
|
wolffd@0
|
23 % 1) a root node
|
wolffd@0
|
24 % 2) linked to exactly two observed variables.
|
wolffd@0
|
25 % Latent variables in the projection are represented using a bidirectional graph,
|
wolffd@0
|
26 % and thus remain implicit.
|
wolffd@0
|
27 %
|
wolffd@0
|
28 % See Pearl, "Causality: Models, Reasoning, and Inference", 2000, p52 for more details.
|
wolffd@0
|
29 % Written by Tamar Kushnir, 2000
|
wolffd@0
|
30
|
wolffd@0
|
31 sep = cell(n,n);
|
wolffd@0
|
32 ord = 0;
|
wolffd@0
|
33 done = 0;
|
wolffd@0
|
34 G = ones(n,n);
|
wolffd@0
|
35 G = setdiag(G,0);
|
wolffd@0
|
36 while ~done
|
wolffd@0
|
37 done = 1;
|
wolffd@0
|
38 [X,Y] = find(G);
|
wolffd@0
|
39 for i=1:length(X)
|
wolffd@0
|
40 x = X(i); y = Y(i);
|
wolffd@0
|
41 nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]);
|
wolffd@0
|
42 if length(nbrs) >= ord & G(x,y) ~= 0
|
wolffd@0
|
43 done = 0;
|
wolffd@0
|
44 SS = subsets(nbrs, ord, ord); % all subsets of size ord
|
wolffd@0
|
45 for si=1:length(SS)
|
wolffd@0
|
46 S = SS{si};
|
wolffd@0
|
47 if feval(cond_indep, x, y, S, varargin{:})
|
wolffd@0
|
48 G(x,y) = 0;
|
wolffd@0
|
49 G(y,x) = 0;
|
wolffd@0
|
50 sep{x,y} = myunion(sep{x,y}, S);
|
wolffd@0
|
51 sep{y,x} = myunion(sep{y,x}, S);
|
wolffd@0
|
52 break; % no need to check any more subsets
|
wolffd@0
|
53 end
|
wolffd@0
|
54 end
|
wolffd@0
|
55 end
|
wolffd@0
|
56 end
|
wolffd@0
|
57 ord = ord + 1;
|
wolffd@0
|
58 end
|
wolffd@0
|
59
|
wolffd@0
|
60 % Create the minimal pattern,
|
wolffd@0
|
61 % i.e., the only directed edges are V structures.
|
wolffd@0
|
62 pdag = G;
|
wolffd@0
|
63 [X, Y] = find(G);
|
wolffd@0
|
64 % We want to generate all unique triples x,y,z
|
wolffd@0
|
65 % where y is a common neighbor to x and z
|
wolffd@0
|
66 for i=1:length(X)
|
wolffd@0
|
67 x = X(i);
|
wolffd@0
|
68 y = Y(i);
|
wolffd@0
|
69 Z = find(G(y,:));
|
wolffd@0
|
70 Z = mysetdiff(Z, x);
|
wolffd@0
|
71 for z=Z(:)'
|
wolffd@0
|
72 if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x})
|
wolffd@0
|
73 pdag(x,y) = -1; pdag(y,x) = 0;
|
wolffd@0
|
74 pdag(z,y) = -1; pdag(y,z) = 0;
|
wolffd@0
|
75 end
|
wolffd@0
|
76 end
|
wolffd@0
|
77 end
|
wolffd@0
|
78
|
wolffd@0
|
79 % Convert the minimal pattern to a complete one using the following rules:
|
wolffd@0
|
80 % Rule 1:
|
wolffd@0
|
81 % if a and b are non-adjacent nodes with a common neighbor c,
|
wolffd@0
|
82 % if a->c and not b->c then c-*>b (marked arrow).
|
wolffd@0
|
83 % Rule 2:
|
wolffd@0
|
84 % if a and b are adjacent and there is a directed path (marked links) from a to b
|
wolffd@0
|
85 % then a->b (add arrowhead).
|
wolffd@0
|
86 %Pearl (2000)
|
wolffd@0
|
87
|
wolffd@0
|
88 arrowin = [-1 -2 2];
|
wolffd@0
|
89 old_pdag = zeros(n);
|
wolffd@0
|
90 iter = 0;
|
wolffd@0
|
91 while ~isequal(pdag, old_pdag)
|
wolffd@0
|
92 iter = iter + 1;
|
wolffd@0
|
93 old_pdag = pdag;
|
wolffd@0
|
94 % rule 1
|
wolffd@0
|
95 [X, Y] = find(pdag);
|
wolffd@0
|
96 for i=1:length(X)
|
wolffd@0
|
97 x = X(i);
|
wolffd@0
|
98 y = Y(i);
|
wolffd@0
|
99 Z = find(pdag(y,:));
|
wolffd@0
|
100 Z = mysetdiff(Z, x);
|
wolffd@0
|
101 for z=Z(:)'
|
wolffd@0
|
102 if G(x,z)==0 & ismember(pdag(x,y),arrowin) & ~ismember(pdag(z,y),arrowin)
|
wolffd@0
|
103 pdag(y,z) = -2; pdag(z,y) = 0;
|
wolffd@0
|
104 end
|
wolffd@0
|
105 end
|
wolffd@0
|
106 end
|
wolffd@0
|
107 % rule 2
|
wolffd@0
|
108 [X, Y] = find(G);
|
wolffd@0
|
109 %check all adjacent nodes because if pdag(x,y) = -1
|
wolffd@0
|
110 %and pdag(y,x) = 0 there could still be an bidirected edge between x & y.
|
wolffd@0
|
111 for i=1:length(X)
|
wolffd@0
|
112 x = X(i);
|
wolffd@0
|
113 y = Y(i);
|
wolffd@0
|
114 if ~ismember(pdag(x,y), arrowin) %x->y doesn't exist yet
|
wolffd@0
|
115 %find marked path from x to y
|
wolffd@0
|
116 add_arrow = marked_path(x,y,pdag);
|
wolffd@0
|
117 if add_arrow
|
wolffd@0
|
118 if pdag(y,x)==-1 %bidirected edge
|
wolffd@0
|
119 pdag(x,y) = 2; pdag(y,x) = 2;
|
wolffd@0
|
120 else
|
wolffd@0
|
121 pdag(x,y) = -1;pdag(y,x) = 0;
|
wolffd@0
|
122 end
|
wolffd@0
|
123 end
|
wolffd@0
|
124 end
|
wolffd@0
|
125 end
|
wolffd@0
|
126 end
|
wolffd@0
|
127
|
wolffd@0
|
128
|
wolffd@0
|
129 %%%%%%%%%%%%%
|
wolffd@0
|
130
|
wolffd@0
|
131 function t = marked_path(x,y,L)
|
wolffd@0
|
132 % MARKED_PATH is a boolean function which returns 1 if a marked path
|
wolffd@0
|
133 % between nodes x and y exists in the partially directed latent structure L.
|
wolffd@0
|
134 %
|
wolffd@0
|
135 % t = marked_path(x,y,L)
|
wolffd@0
|
136 %
|
wolffd@0
|
137 % x and y are the starting and ending nodes in the path, respectively.
|
wolffd@0
|
138 % L is a latent structure (partially directed graph with possible latent variables).
|
wolffd@0
|
139 %
|
wolffd@0
|
140 % Rule 2 of IC* algorithm (see Pearl, 2000)
|
wolffd@0
|
141
|
wolffd@0
|
142 t=0;
|
wolffd@0
|
143
|
wolffd@0
|
144 %find set of marked links from x
|
wolffd@0
|
145 marked = find(L(x,:)==-2);
|
wolffd@0
|
146 if ismember(y,marked)
|
wolffd@0
|
147 t=1; %marked path found
|
wolffd@0
|
148 else
|
wolffd@0
|
149 for m=marked(:)'
|
wolffd@0
|
150 t = marked_path(m,y,L);
|
wolffd@0
|
151 if t==1
|
wolffd@0
|
152 break; %stop when marked path found
|
wolffd@0
|
153 end
|
wolffd@0
|
154 end
|
wolffd@0
|
155 end
|