Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/general/convert_dbn_CPDs_to_tables.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 CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence)
wolffd@0 2 % CONVERT_DBN_CPDS_TO_TABLES Convert CPDs of (possibly instantiated) DBN nodes to tables
wolffd@0 3 % CPDpot = convert_dbn_CPDs_to_tables(bnet, evidence)
wolffd@0 4 %
wolffd@0 5 % CPDpot{n,t} is a table containing P(n,t|pa(n,t), ev)
wolffd@0 6 % All hidden nodes are assumed to be discrete.
wolffd@0 7 % We assume the observed nodes are the same in every slice.
wolffd@0 8 %
wolffd@0 9 % Evaluating the conditional likelihood of long evidence sequences can be very slow,
wolffd@0 10 % so we take pains to vectorize where possible.
wolffd@0 11
wolffd@0 12 [ss T] = size(evidence);
wolffd@0 13 %obs_bitv = ~isemptycell(evidence(:));
wolffd@0 14 obs_bitv = zeros(1, 2*ss);
wolffd@0 15 obs_bitv(bnet.observed) = 1;
wolffd@0 16 obs_bitv(bnet.observed+ss) = 1;
wolffd@0 17
wolffd@0 18 ns = bnet.node_sizes(:);
wolffd@0 19 CPDpot = cell(ss,T);
wolffd@0 20
wolffd@0 21 for n=1:ss
wolffd@0 22 % slice 1
wolffd@0 23 t = 1;
wolffd@0 24 ps = parents(bnet.dag, n);
wolffd@0 25 e = bnet.equiv_class(n, 1);
wolffd@0 26 if ~any(obs_bitv(ps))
wolffd@0 27 CPDpot{n,t} = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, evidence{n,t});
wolffd@0 28 else
wolffd@0 29 CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps n], evidence(:,1));
wolffd@0 30 end
wolffd@0 31
wolffd@0 32 % special cases: c=child, p=parents, d=discrete, h=hidden, 1sl=1slice
wolffd@0 33 % if c=h=1 then c=d=1, since hidden nodes must be discrete
wolffd@0 34 % c=h c=d p=h p=d 1sl method
wolffd@0 35 % ---------------------------
wolffd@0 36 % 1 1 1 1 - replicate CPT
wolffd@0 37 % - 1 - 1 - evaluate CPT on evidence *
wolffd@0 38 % 0 1 1 1 1 dhmm
wolffd@0 39 % 0 0 1 1 1 ghmm
wolffd@0 40 % other loop
wolffd@0 41 %
wolffd@0 42 % * = any subset of the domain may be observed
wolffd@0 43
wolffd@0 44 % Example where all of the special cases occur - a hierarchical HMM
wolffd@0 45 % where the top layer (G) and leaves (Y) are observed and
wolffd@0 46 % all nodes are discrete except Y.
wolffd@0 47 % (O turns on if Y is an outlier)
wolffd@0 48
wolffd@0 49 % G ---------> G
wolffd@0 50 % | |
wolffd@0 51 % v v
wolffd@0 52 % S --------> S
wolffd@0 53 % | |
wolffd@0 54 % v v
wolffd@0 55 % Y Y
wolffd@0 56 % ^ ^
wolffd@0 57 % | |
wolffd@0 58 % O O
wolffd@0 59
wolffd@0 60 % Evaluating P(yt|St,Ot) is the ghmm case
wolffd@0 61 % Evaluating P(St|S(t-1),gt) is the eval CPT case
wolffd@0 62 % Evaluating P(gt|g(t-1) is the eval CPT case (hdom = [])
wolffd@0 63 % Evaluating P(Ot) is the replicated CPT case
wolffd@0 64
wolffd@0 65 % Cts parents (e.g., inputs) would require an additional special case for speed
wolffd@0 66
wolffd@0 67
wolffd@0 68 % slices 2..T
wolffd@0 69 [ss T] = size(evidence);
wolffd@0 70 self = n+ss;
wolffd@0 71 ps = parents(bnet.dag, self);
wolffd@0 72 e = bnet.equiv_class(n, 2);
wolffd@0 73
wolffd@0 74 if 1
wolffd@0 75 debug = 0;
wolffd@0 76 hidden_child = ~obs_bitv(n);
wolffd@0 77 discrete_child = myismember(n, bnet.dnodes);
wolffd@0 78 hidden_ps = all(~obs_bitv(ps));
wolffd@0 79 discrete_ps = mysubset(ps, bnet.dnodes);
wolffd@0 80 parents_in_same_slice = all(ps > ss);
wolffd@0 81
wolffd@0 82 if hidden_child & discrete_child & hidden_ps & discrete_ps
wolffd@0 83 CPDpot = helper_repl(bnet, evidence, n, CPDpot, obs_bitv, debug);
wolffd@0 84 elseif discrete_child & discrete_ps
wolffd@0 85 CPDpot = helper_eval(bnet, evidence, n, CPDpot, obs_bitv, debug);
wolffd@0 86 elseif discrete_child & hidden_ps & discrete_ps & parents_in_same_slice
wolffd@0 87 CPDpot = helper_dhmm(bnet, evidence, n, CPDpot, obs_bitv, debug);
wolffd@0 88 elseif ~discrete_child & hidden_ps & discrete_ps & parents_in_same_slice
wolffd@0 89 CPDpot = helper_ghmm(bnet, evidence, n, CPDpot, obs_bitv, debug);
wolffd@0 90 else
wolffd@0 91 if debug, fprintf('node %d, slow\n', n); end
wolffd@0 92 for t=2:T
wolffd@0 93 CPDpot{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
wolffd@0 94 end
wolffd@0 95 end
wolffd@0 96 end
wolffd@0 97
wolffd@0 98 if 0
wolffd@0 99 for t=2:T
wolffd@0 100 CPDpot2{n,t} = convert_to_table(bnet.CPD{e}, [ps self], evidence(:,t-1:t));
wolffd@0 101 if ~approxeq(CPDpot{n,t}, CPDpot2{n,t})
wolffd@0 102 fprintf('CPDpot n=%d, t=%d\n',n,t);
wolffd@0 103 keyboard
wolffd@0 104 end
wolffd@0 105 end
wolffd@0 106 end
wolffd@0 107
wolffd@0 108
wolffd@0 109 end
wolffd@0 110
wolffd@0 111
wolffd@0 112
wolffd@0 113
wolffd@0 114 %%%%%%%
wolffd@0 115 function CPDpot = helper_repl(bnet, evidence, n, CPDpot, obs_bitv, debug)
wolffd@0 116
wolffd@0 117 [ss T] = size(evidence);
wolffd@0 118 if debug, fprintf('node %d, repl\n', n); end
wolffd@0 119 e = bnet.equiv_class(n, 2);
wolffd@0 120 CPT = convert_CPD_to_table_hidden_ps(bnet.CPD{e}, []);
wolffd@0 121 CPDpot(n,2:T) = num2cell(repmat(CPT, [1 1 T-1]), [1 2]);
wolffd@0 122
wolffd@0 123
wolffd@0 124
wolffd@0 125 %%%%%%%
wolffd@0 126 function CPDpot = helper_eval(bnet, evidence, n, CPDpot, obs_bitv, debug)
wolffd@0 127
wolffd@0 128 [ss T] = size(evidence);
wolffd@0 129 self = n+ss;
wolffd@0 130 ps = parents(bnet.dag, self);
wolffd@0 131 e = bnet.equiv_class(n, 2);
wolffd@0 132 ns = bnet.node_sizes(:);
wolffd@0 133 % Example: given CPT(p1, p2, p3, p4, c), where p1,p3 are observed
wolffd@0 134 % we create CPT([p2 p4 c], [p1 p3]).
wolffd@0 135 % We then convert all observed p1,p3 into indices ndx
wolffd@0 136 % and return CPT(:, ndx)
wolffd@0 137 CPT = CPD_to_CPT(bnet.CPD{e});
wolffd@0 138 domain = [ps self];
wolffd@0 139 % if dom is [3 7 8] and 3,8 are observed, odom_rel = [1 3], hdom_rel = 2,
wolffd@0 140 % odom = [3 8], hdom = 7
wolffd@0 141 odom_rel = find(obs_bitv(domain));
wolffd@0 142 hdom_rel = find(~obs_bitv(domain));
wolffd@0 143 odom = domain(odom_rel);
wolffd@0 144 hdom = domain(hdom_rel);
wolffd@0 145 if isempty(hdom)
wolffd@0 146 CPT = CPT(:);
wolffd@0 147 else
wolffd@0 148 CPT = permute(CPT, [hdom_rel odom_rel]);
wolffd@0 149 CPT = reshape(CPT, prod(ns(hdom)), prod(ns(odom)));
wolffd@0 150 end
wolffd@0 151 parents_in_same_slice = all(ps > ss);
wolffd@0 152 if parents_in_same_slice
wolffd@0 153 if debug, fprintf('node %d eval 1 slice\n', n); end
wolffd@0 154 data = cell2num(evidence(odom-ss,2:T)); %data(i,t) = val of i'th obs parent at t+1
wolffd@0 155 else
wolffd@0 156 if debug, fprintf('node %d eval 2 slice\n', n); end
wolffd@0 157 % there's probably a way of vectorizing this...
wolffd@0 158 data = zeros(length(odom), T-1);
wolffd@0 159 for t=2:T
wolffd@0 160 ev = evidence(:,t-1:t);
wolffd@0 161 ev = ev(:);
wolffd@0 162 ev2 = ev(odom);
wolffd@0 163 data(:,t-1) = cat(1, ev2{:});
wolffd@0 164 %data(:,t-1) = cell2num(ev2);
wolffd@0 165 end
wolffd@0 166 end
wolffd@0 167 ndx = subv2ind(ns(odom), data'); % ndx(t) encodes data(:,t)
wolffd@0 168 if isempty(hdom)
wolffd@0 169 CPDpot(n,2:T) = num2cell(CPT(ndx)); % a cell array of floats
wolffd@0 170 else
wolffd@0 171 CPDpot(n,2:T) = num2cell(CPT(:, ndx), 1); % a cell array of column vectors
wolffd@0 172 end
wolffd@0 173
wolffd@0 174 %%%%%%%
wolffd@0 175 function CPDpot = helper_dhmm(bnet, evidence, n, CPDpot, obs_bitv, debug)
wolffd@0 176
wolffd@0 177 if debug, fprintf('node %d, dhmm\n', n); end
wolffd@0 178 [ss T] = size(evidence);
wolffd@0 179 self = n+ss;
wolffd@0 180 ps = parents(bnet.dag, self);
wolffd@0 181 e = bnet.equiv_class(n, 2);
wolffd@0 182 ns = bnet.node_sizes(:);
wolffd@0 183 CPT = CPD_to_CPT(bnet.CPD{e});
wolffd@0 184 CPT = reshape(CPT, [prod(ns(ps)) ns(self)]); % what if no parents?
wolffd@0 185 %obslik = mk_dhmm_obs_lik(cell2num(evidence(n,2:T)), CPT);
wolffd@0 186 obslik = eval_pdf_cond_multinomial(cell2num(evidence(n,2:T)), CPT);
wolffd@0 187 CPDpot(n,2:T) = num2cell(obslik, 1);
wolffd@0 188
wolffd@0 189
wolffd@0 190 %%%%%%%
wolffd@0 191 function CPDpot = helper_ghmm(bnet, evidence, n, CPDpot, obs_bitv, debug)
wolffd@0 192
wolffd@0 193 if debug, fprintf('node %d, ghmm\n', n); end
wolffd@0 194 [ss T] = size(evidence);
wolffd@0 195 e = bnet.equiv_class(n, 2);
wolffd@0 196 S = struct(bnet.CPD{e});
wolffd@0 197 ev2 = cell2num(evidence(n,2:T));
wolffd@0 198 %obslik = mk_ghmm_obs_lik(ev2, S.mean, S.cov);
wolffd@0 199 obslik = eval_pdf_cond_gauss(ev2, S.mean, S.cov);
wolffd@0 200 CPDpot(n,2:T) = num2cell(obslik, 1);
wolffd@0 201