Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/learning/learn_struct_pdag_pc.m @ 0:e9a9cd732c1e tip

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
first hg version after svn
author wolffd
date Tue, 10 Feb 2015 15:05:51 +0000
parents
children
rev   line source
wolffd@0 1 function [pdag, G] = learn_struct_pdag_pc(cond_indep, n, k, varargin)
wolffd@0 2 % LEARN_STRUCT_PDAG_PC Learn a partially oriented DAG (pattern) using the PC algorithm
wolffd@0 3 % P = learn_struct_pdag_pc(cond_indep, n, k, ...)
wolffd@0 4 %
wolffd@0 5 % n is the number of nodes.
wolffd@0 6 % k is an optional upper bound on the fan-in (default: n)
wolffd@0 7 % cond_indep is a boolean function that will be called as follows:
wolffd@0 8 % feval(cond_indep, x, y, S, ...)
wolffd@0 9 % where x and y are nodes, and S is a set of nodes (positive integers),
wolffd@0 10 % and ... are any optional parameters passed to this function.
wolffd@0 11 %
wolffd@0 12 % The output P is an adjacency matrix, in which
wolffd@0 13 % P(i,j) = -1 if there is an i->j edge.
wolffd@0 14 % P(i,j) = P(j,i) = 1 if there is an undirected edge i <-> j
wolffd@0 15 %
wolffd@0 16 % The PC algorithm does structure learning assuming all variables are observed.
wolffd@0 17 % See Spirtes, Glymour and Scheines, "Causation, Prediction and Search", 1993, p117.
wolffd@0 18 % This algorithm may take O(n^k) time if there are n variables and k is the max fan-in,
wolffd@0 19 % but this is quicker than the Verma-Pearl IC algorithm, which is always O(n^n).
wolffd@0 20
wolffd@0 21
wolffd@0 22 sep = cell(n,n);
wolffd@0 23 ord = 0;
wolffd@0 24 done = 0;
wolffd@0 25 G = ones(n,n);
wolffd@0 26 G=setdiag(G,0);
wolffd@0 27 while ~done
wolffd@0 28 done = 1;
wolffd@0 29 [X,Y] = find(G);
wolffd@0 30 for i=1:length(X)
wolffd@0 31 x = X(i); y = Y(i);
wolffd@0 32 %nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]);
wolffd@0 33 nbrs = mysetdiff(neighbors(G, y), x); % bug fix by Raanan Yehezkel <raanany@ee.bgu.ac.il> 6/27/04
wolffd@0 34 if length(nbrs) >= ord & G(x,y) ~= 0
wolffd@0 35 done = 0;
wolffd@0 36 %SS = subsets(nbrs, ord, ord); % all subsets of size ord
wolffd@0 37 SS = subsets1(nbrs, ord);
wolffd@0 38 for si=1:length(SS)
wolffd@0 39 S = SS{si};
wolffd@0 40 if feval(cond_indep, x, y, S, varargin{:})
wolffd@0 41 %if isempty(S)
wolffd@0 42 % fprintf('%d indep of %d ', x, y);
wolffd@0 43 %else
wolffd@0 44 % fprintf('%d indep of %d given ', x, y); fprintf('%d ', S);
wolffd@0 45 %end
wolffd@0 46 %fprintf('\n');
wolffd@0 47
wolffd@0 48 % diagnostic
wolffd@0 49 %[CI, r] = cond_indep_fisher_z(x, y, S, varargin{:});
wolffd@0 50 %fprintf(': r = %6.4f\n', r);
wolffd@0 51
wolffd@0 52 G(x,y) = 0;
wolffd@0 53 G(y,x) = 0;
wolffd@0 54 sep{x,y} = myunion(sep{x,y}, S);
wolffd@0 55 sep{y,x} = myunion(sep{y,x}, S);
wolffd@0 56 break; % no need to check any more subsets
wolffd@0 57 end
wolffd@0 58 end
wolffd@0 59 end
wolffd@0 60 end
wolffd@0 61 ord = ord + 1;
wolffd@0 62 end
wolffd@0 63
wolffd@0 64
wolffd@0 65 % Create the minimal pattern,
wolffd@0 66 % i.e., the only directed edges are V structures.
wolffd@0 67 pdag = G;
wolffd@0 68 [X, Y] = find(G);
wolffd@0 69 % We want to generate all unique triples x,y,z
wolffd@0 70 % This code generates x,y,z and z,y,x.
wolffd@0 71 for i=1:length(X)
wolffd@0 72 x = X(i);
wolffd@0 73 y = Y(i);
wolffd@0 74 Z = find(G(y,:));
wolffd@0 75 Z = mysetdiff(Z, x);
wolffd@0 76 for z=Z(:)'
wolffd@0 77 if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x})
wolffd@0 78 %fprintf('%d -> %d <- %d\n', x, y, z);
wolffd@0 79 pdag(x,y) = -1; pdag(y,x) = 0;
wolffd@0 80 pdag(z,y) = -1; pdag(y,z) = 0;
wolffd@0 81 end
wolffd@0 82 end
wolffd@0 83 end
wolffd@0 84
wolffd@0 85 % Convert the minimal pattern to a complete one,
wolffd@0 86 % i.e., every directed edge in P is compelled
wolffd@0 87 % (must be directed in all Markov equivalent models),
wolffd@0 88 % and every undirected edge in P is reversible.
wolffd@0 89 % We use the rules of Pearl (2000) p51 (derived in Meek (1995))
wolffd@0 90
wolffd@0 91 old_pdag = zeros(n);
wolffd@0 92 iter = 0;
wolffd@0 93 while ~isequal(pdag, old_pdag)
wolffd@0 94 iter = iter + 1;
wolffd@0 95 old_pdag = pdag;
wolffd@0 96 % rule 1
wolffd@0 97 [A,B] = find(pdag==-1); % a -> b
wolffd@0 98 for i=1:length(A)
wolffd@0 99 a = A(i); b = B(i);
wolffd@0 100 C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a
wolffd@0 101 if ~isempty(C)
wolffd@0 102 pdag(b,C) = -1; pdag(C,b) = 0;
wolffd@0 103 %fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C);
wolffd@0 104 end
wolffd@0 105 end
wolffd@0 106 % rule 2
wolffd@0 107 [A,B] = find(pdag==1); % unoriented a-b edge
wolffd@0 108 for i=1:length(A)
wolffd@0 109 a = A(i); b = B(i);
wolffd@0 110 if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' );
wolffd@0 111 pdag(a,b) = -1; pdag(b,a) = 0;
wolffd@0 112 %fprintf('rule 2: %d -> %d\n', a, b);
wolffd@0 113 end
wolffd@0 114 end
wolffd@0 115 % rule 3
wolffd@0 116 [A,B] = find(pdag==1); % a-b
wolffd@0 117 for i=1:length(A)
wolffd@0 118 a = A(i); b = B(i);
wolffd@0 119 C = find( (pdag(a,:)==1) & (pdag(:,b)==-1)' );
wolffd@0 120 % C contains nodes c s.t. a-c->ba
wolffd@0 121 G2 = setdiag(G(C, C), 1);
wolffd@0 122 if any(G2(:)==0) % there are 2 different non adjacent elements of C
wolffd@0 123 pdag(a,b) = -1; pdag(b,a) = 0;
wolffd@0 124 %fprintf('rule 3: %d -> %d\n', a, b);
wolffd@0 125 end
wolffd@0 126 end
wolffd@0 127 end
wolffd@0 128
wolffd@0 129