Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/learning/learn_struct_pdag_pc_constrain.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] = dn_learn_struct_pdag_pc_constrain(adj, 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 % Pdag = learn_struct_pdag_pc_constrain(adj, cond_indep, n, k, ...)
wolffd@0 4 %
wolffd@0 5 % adj = adjacency matrix learned from dependency network P(i,j) = 1 => i--j; 0 => i j
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 %Output
wolffd@0 14 % pdag Partially directed graph
wolffd@0 15 % G Resulting adjacency graph prior to setting direction arrows
wolffd@0 16 %
wolffd@0 17 % The output P is an adjacency matrix, in which
wolffd@0 18 % P(i,j) = -1 if there is an i->j edge.
wolffd@0 19 % P(i,j) = P(j,i) = 1 if there is an undirected edge i <-> j
wolffd@0 20 %
wolffd@0 21 % The PC algorithm does structure learning assuming all variables are observed.
wolffd@0 22 % See Spirtes, Glymour and Scheines, "Causation, Prediction and Search", 1993, p117.
wolffd@0 23 % This algorithm may take O(n^k) time if there are n variables and k is the max fan-in,
wolffd@0 24 % but this is quicker than the Verma-Pearl IC algorithm, which is always O(n^n).
wolffd@0 25 %
wolffd@0 26 %% Example
wolffd@0 27 %% Given data in a comma separated, filename starting with the variable labels, then the data in rows.
wolffd@0 28 %% filename test.txt consists of:
wolffd@0 29 %%
wolffd@0 30 %% Earthquake,Burglar,Radio,Alarm,Call
wolffd@0 31 %% 1,2,2,2,1
wolffd@0 32 %% 1,1,2,1,2
wolffd@0 33 %% . . .
wolffd@0 34 %[CovMatrix, obs, varfields] = CovMat('test.txt',5);
wolffd@0 35 %
wolffd@0 36 %dn = zeros(5,5);
wolffd@0 37 %dn(1,2) = 1; % This was the known Markov blanket of the system that generated test.txt
wolffd@0 38 %dn(2,1) = 1;
wolffd@0 39 %dn(2,4) = 1;
wolffd@0 40 %dn(4,2) = 1;
wolffd@0 41 %dn(1,3) = 1;
wolffd@0 42 %dn(3,1) = 1;
wolffd@0 43 %dn(1,4) = 1;
wolffd@0 44 %dn(4,1) = 1;
wolffd@0 45 %dn(4,5) = 1;
wolffd@0 46 %dn(5,4) = 1;
wolffd@0 47 %dn(3,5) = 1; %loop r->c
wolffd@0 48 %dn(5,3) = 1; %loop c-r
wolffd@0 49 %dn(3,4) = 1;
wolffd@0 50 %dn(4,3) = 1;
wolffd@0 51 %
wolffd@0 52 %max_fan_in = 4;
wolffd@0 53 %alpha = 0.05;
wolffd@0 54 %
wolffd@0 55 %[pdag G] = learn_struct_pdag_pc_constrain(dn,'cond_indep_fisher_z', 5, max_fan_in, CovMatrix, obs, alpha);
wolffd@0 56 %%
wolffd@0 57 %%
wolffd@0 58 %% Gary Bradski, 7/2002 Modified this to take an adjacency matrix from a dependency network.
wolffd@0 59
wolffd@0 60
wolffd@0 61 sep = cell(n,n);
wolffd@0 62 ord = 0;
wolffd@0 63 done = 0;
wolffd@0 64 G = ones(n,n);
wolffd@0 65 G=setdiag(G,0);
wolffd@0 66
wolffd@0 67 while ~done
wolffd@0 68 done = 1;
wolffd@0 69 [X,Y] = find(G);
wolffd@0 70 for i=1:length(X)
wolffd@0 71 x = X(i); y = Y(i);
wolffd@0 72 % nbrs = mysetdiff(myunion(neighbors(G, x), neighbors(G,y)), [x y]);%parents, children, but not self
wolffd@0 73 nbrs = mysetdiff(myunion(neighbors(adj, x), neighbors(adj,y)), [x y]);%parents, children, but not self
wolffd@0 74
wolffd@0 75 if length(nbrs) >= ord & G(x,y) ~= 0
wolffd@0 76 done = 0;
wolffd@0 77 SS = subsets(nbrs, ord, ord); % all subsets of size ord
wolffd@0 78 for si=1:length(SS)
wolffd@0 79 S = SS{si};
wolffd@0 80 %if (feval(dsep,x,y,S,adj)) | (feval(cond_indep, x, y, S, varargin{:}))
wolffd@0 81 if feval(cond_indep, x, y, S, varargin{:})
wolffd@0 82 %if isempty(S)
wolffd@0 83 % fprintf('%d indep of %d ', x, y);
wolffd@0 84 %else
wolffd@0 85 % fprintf('%d indep of %d given ', x, y); fprintf('%d ', S);
wolffd@0 86 %end
wolffd@0 87 %fprintf('\n');
wolffd@0 88
wolffd@0 89 % diagnostic
wolffd@0 90 %[CI, r] = cond_indep_fisher_z(x, y, S, varargin{:});
wolffd@0 91 %fprintf(': r = %6.4f\n', r);
wolffd@0 92
wolffd@0 93 G(x,y) = 0;
wolffd@0 94 G(y,x) = 0;
wolffd@0 95 adj(x,y) = 0; %make sure found cond. independencies are marked out
wolffd@0 96 adj(y,x) = 0;
wolffd@0 97 sep{x,y} = myunion(sep{x,y}, S);
wolffd@0 98 sep{y,x} = myunion(sep{y,x}, S);
wolffd@0 99 break; % no need to check any more subsets
wolffd@0 100 end
wolffd@0 101 end
wolffd@0 102 end
wolffd@0 103 end
wolffd@0 104 ord = ord + 1;
wolffd@0 105 end
wolffd@0 106
wolffd@0 107
wolffd@0 108
wolffd@0 109
wolffd@0 110 % Create the minimal pattern,
wolffd@0 111 % i.e., the only directed edges are V structures.
wolffd@0 112
wolffd@0 113 pdag = G;
wolffd@0 114 [X, Y] = find(G);
wolffd@0 115 % We want to generate all unique triples x,y,z
wolffd@0 116 % This code generates x,y,z and z,y,x.
wolffd@0 117 for i=1:length(X)
wolffd@0 118 x = X(i);
wolffd@0 119 y = Y(i);
wolffd@0 120 Z = find(G(y,:));
wolffd@0 121 Z = mysetdiff(Z, x);
wolffd@0 122 for z=Z(:)'
wolffd@0 123 if G(x,z)==0 & ~ismember(y, sep{x,z}) & ~ismember(y, sep{z,x})
wolffd@0 124 %fprintf('%d -> %d <- %d\n', x, y, z);
wolffd@0 125 pdag(x,y) = -1; pdag(y,x) = 0;
wolffd@0 126 pdag(z,y) = -1; pdag(y,z) = 0;
wolffd@0 127 end
wolffd@0 128 end
wolffd@0 129 end
wolffd@0 130
wolffd@0 131 % Convert the minimal pattern to a complete one,
wolffd@0 132 % i.e., every directed edge in P is compelled
wolffd@0 133 % (must be directed in all Markov equivalent models),
wolffd@0 134 % and every undirected edge in P is reversible.
wolffd@0 135 % We use the rules of Pearl (2000) p51 (derived in Meek (1995))
wolffd@0 136
wolffd@0 137 old_pdag = zeros(n);
wolffd@0 138 iter = 0;
wolffd@0 139 while ~isequal(pdag, old_pdag)
wolffd@0 140 iter = iter + 1;
wolffd@0 141 old_pdag = pdag;
wolffd@0 142 % rule 1
wolffd@0 143 [A,B] = find(pdag==-1); % a -> b
wolffd@0 144 for i=1:length(A)
wolffd@0 145 a = A(i); b = B(i);
wolffd@0 146 C = find(pdag(b,:)==1 & G(a,:)==0); % all nodes adj to b but not a
wolffd@0 147 if ~isempty(C)
wolffd@0 148 pdag(b,C) = -1; pdag(C,b) = 0;
wolffd@0 149 %fprintf('rule 1: a=%d->b=%d and b=%d-c=%d implies %d->%d\n', a, b, b, C, b, C);
wolffd@0 150 end
wolffd@0 151 end
wolffd@0 152 % rule 2
wolffd@0 153 [A,B] = find(pdag==1); % unoriented a-b edge
wolffd@0 154 for i=1:length(A)
wolffd@0 155 a = A(i); b = B(i);
wolffd@0 156 if any( (pdag(a,:)==-1) & (pdag(:,b)==-1)' );
wolffd@0 157 pdag(a,b) = -1; pdag(b,a) = 0;
wolffd@0 158 %fprintf('rule 2: %d -> %d\n', a, b);
wolffd@0 159 end
wolffd@0 160 end
wolffd@0 161 % rule 3
wolffd@0 162 [A,B] = find(pdag==1); % a-b
wolffd@0 163 for i=1:length(A)
wolffd@0 164 a = A(i); b = B(i);
wolffd@0 165 C = find( (G(a,:)==1) & (pdag(:,b)==-1)' );
wolffd@0 166 % C contains nodes c s.t. a-c->ba
wolffd@0 167 G2 = setdiag(G(C, C), 1);
wolffd@0 168 if any(G2(:)==0) % there are 2 different non adjacent elements of C
wolffd@0 169 pdag(a,b) = -1; pdag(b,a) = 0;
wolffd@0 170 %fprintf('rule 3: %d -> %d\n', a, b);
wolffd@0 171 end
wolffd@0 172 end
wolffd@0 173 end
wolffd@0 174
wolffd@0 175
wolffd@0 176
wolffd@0 177