Mercurial > hg > camir-aes2014
view 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 |
line wrap: on
line source
function [pdag, G] = learn_struct_pdag_ic_star(cond_indep, n, k, varargin) % LEARN_STRUCT_PDAG_IC_STAR Learn a partially oriented DAG (pattern) with latent % variables using the IC* algorithm % P = learn_struct_pdag_ic_star(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 either a latent variable L such that i <-L-> j % OR there is a directed edge from i->j. % P(i,j) = -2 if there is a marked directed i-*>j edge. % P(i,j) = P(j,i) = 1 if there is and undirected edge i--j % P(i,j) = P(j,i) = 2 if there is a latent variable L such that i<-L->j. % % The IC* algorithm learns a latent structure associated with a set of observed % variables. % The latent structure revealed is the projection in which every latent variable is % 1) a root node % 2) linked to exactly two observed variables. % Latent variables in the projection are represented using a bidirectional graph, % and thus remain implicit. % % See Pearl, "Causality: Models, Reasoning, and Inference", 2000, p52 for more details. % Written by Tamar Kushnir, 2000 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]); if length(nbrs) >= ord & G(x,y) ~= 0 done = 0; SS = subsets(nbrs, ord, ord); % all subsets of size ord for si=1:length(SS) S = SS{si}; if feval(cond_indep, x, y, S, varargin{:}) 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 % where y is a common neighbor to x and z 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}) 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 using the following rules: % Rule 1: % if a and b are non-adjacent nodes with a common neighbor c, % if a->c and not b->c then c-*>b (marked arrow). % Rule 2: % if a and b are adjacent and there is a directed path (marked links) from a to b % then a->b (add arrowhead). %Pearl (2000) arrowin = [-1 -2 2]; old_pdag = zeros(n); iter = 0; while ~isequal(pdag, old_pdag) iter = iter + 1; old_pdag = pdag; % rule 1 [X, Y] = find(pdag); for i=1:length(X) x = X(i); y = Y(i); Z = find(pdag(y,:)); Z = mysetdiff(Z, x); for z=Z(:)' if G(x,z)==0 & ismember(pdag(x,y),arrowin) & ~ismember(pdag(z,y),arrowin) pdag(y,z) = -2; pdag(z,y) = 0; end end end % rule 2 [X, Y] = find(G); %check all adjacent nodes because if pdag(x,y) = -1 %and pdag(y,x) = 0 there could still be an bidirected edge between x & y. for i=1:length(X) x = X(i); y = Y(i); if ~ismember(pdag(x,y), arrowin) %x->y doesn't exist yet %find marked path from x to y add_arrow = marked_path(x,y,pdag); if add_arrow if pdag(y,x)==-1 %bidirected edge pdag(x,y) = 2; pdag(y,x) = 2; else pdag(x,y) = -1;pdag(y,x) = 0; end end end end end %%%%%%%%%%%%% function t = marked_path(x,y,L) % MARKED_PATH is a boolean function which returns 1 if a marked path % between nodes x and y exists in the partially directed latent structure L. % % t = marked_path(x,y,L) % % x and y are the starting and ending nodes in the path, respectively. % L is a latent structure (partially directed graph with possible latent variables). % % Rule 2 of IC* algorithm (see Pearl, 2000) t=0; %find set of marked links from x marked = find(L(x,:)==-2); if ismember(y,marked) t=1; %marked path found else for m=marked(:)' t = marked_path(m,y,L); if t==1 break; %stop when marked path found end end end