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