Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/inference/dynamic/@pearl_dbn_inf_engine/Old/filter_evidence_obj_oriented.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 [marginal, msg, loglik] = filter_evidence_old(engine, evidence)
wolffd@0 2 % [marginal, msg, loglik] = filter_evidence(engine, evidence) (pearl_dbn)
wolffd@0 3
wolffd@0 4 [ss T] = size(evidence);
wolffd@0 5 bnet = bnet_from_engine(engine);
wolffd@0 6 bnet2 = dbn_to_bnet(bnet, T);
wolffd@0 7 ns = bnet2.node_sizes;
wolffd@0 8 hnodes = mysetdiff(1:ss, engine.onodes);
wolffd@0 9 hnodes = hnodes(:)';
wolffd@0 10
wolffd@0 11 [engine.parent_index, engine.child_index] = mk_pearl_msg_indices(bnet2);
wolffd@0 12
wolffd@0 13 msg = init_msgs(bnet2.dag, ns, evidence);
wolffd@0 14 msg = init_ev_msgs(engine, evidence, msg);
wolffd@0 15
wolffd@0 16 verbose = 1;
wolffd@0 17 if verbose, fprintf('\nold filtering\n'); end
wolffd@0 18
wolffd@0 19 for t=1:T
wolffd@0 20 % update pi
wolffd@0 21 for i=hnodes
wolffd@0 22 n = i + (t-1)*ss;
wolffd@0 23 ps = parents(bnet2.dag, n);
wolffd@0 24 if t==1
wolffd@0 25 e = bnet.equiv_class(i,1);
wolffd@0 26 else
wolffd@0 27 e = bnet.equiv_class(i,2);
wolffd@0 28 end
wolffd@0 29 msg{n}.pi = compute_pi(bnet.CPD{e}, n, ps, msg);
wolffd@0 30 %if verbose, fprintf('%d computes pi\n', n); disp(msg{n}.pi); end
wolffd@0 31 msg{n}.pi = normalise(msg{n}.pi(:) .* msg{n}.lambda_from_self(:));
wolffd@0 32 if verbose, fprintf('%d recomputes pi\n', n); disp(msg{n}.pi); end
wolffd@0 33 end
wolffd@0 34 % send pi msg to children
wolffd@0 35 for i=hnodes
wolffd@0 36 n = i + (t-1)*ss;
wolffd@0 37 cs = children(bnet2.dag, n);
wolffd@0 38 for c=cs(:)'
wolffd@0 39 j = engine.parent_index{c}(n); % n is c's j'th parent
wolffd@0 40 pi_msg = normalise(compute_pi_msg(n, cs, msg, c, ns));
wolffd@0 41 msg{c}.pi_from_parent{j} = pi_msg;
wolffd@0 42 if verbose, fprintf('%d sends pi to %d\n', n,c); disp(pi_msg); end
wolffd@0 43 end
wolffd@0 44 end
wolffd@0 45 end
wolffd@0 46
wolffd@0 47
wolffd@0 48 marginal = cell(ss,T);
wolffd@0 49 lik = zeros(1,ss*T);
wolffd@0 50 for t=1:T
wolffd@0 51 for i=1:ss
wolffd@0 52 n = i + (t-1)*ss;
wolffd@0 53 %[bel, lik(n)] = normalise(msg{n}.pi .* msg{n}.lambda);
wolffd@0 54 [bel, lik(n)] = normalise(msg{n}.pi);
wolffd@0 55 marginal{i,t} = bel;
wolffd@0 56 end
wolffd@0 57 end
wolffd@0 58
wolffd@0 59 loglik = sum(log(lik));
wolffd@0 60
wolffd@0 61
wolffd@0 62
wolffd@0 63 %%%%%%%
wolffd@0 64
wolffd@0 65 function lambda = compute_lambda(n, cs, msg, ns)
wolffd@0 66 % Pearl p183 eq 4.50
wolffd@0 67 lambda = prod_lambda_msgs(n, cs, msg, ns);
wolffd@0 68
wolffd@0 69 %%%%%%%
wolffd@0 70
wolffd@0 71 function pi_msg = compute_pi_msg(n, cs, msg, c, ns)
wolffd@0 72 % Pearl p183 eq 4.53 and 4.51
wolffd@0 73 pi_msg = msg{n}.pi .* prod_lambda_msgs(n, cs, msg, ns, c);
wolffd@0 74
wolffd@0 75 %%%%%%%%%
wolffd@0 76
wolffd@0 77 function lam = prod_lambda_msgs(n, cs, msg, ns, except)
wolffd@0 78
wolffd@0 79 if nargin < 5, except = -1; end
wolffd@0 80
wolffd@0 81 %lam = msg{n}.lambda_from_self(:);
wolffd@0 82 lam = ones(ns(n), 1);
wolffd@0 83 for i=1:length(cs)
wolffd@0 84 c = cs(i);
wolffd@0 85 if c ~= except
wolffd@0 86 lam = lam .* msg{n}.lambda_from_child{i};
wolffd@0 87 end
wolffd@0 88 end
wolffd@0 89
wolffd@0 90
wolffd@0 91 %%%%%%%%%%%
wolffd@0 92
wolffd@0 93 function msg = init_msgs(dag, ns, evidence)
wolffd@0 94 % INIT_MSGS Initialize the lambda/pi message and state vectors (pearl_dbn)
wolffd@0 95 % msg = init_msgs(dag, ns, evidence)
wolffd@0 96 %
wolffd@0 97 % We assume all the hidden nodes are discrete.
wolffd@0 98
wolffd@0 99 N = length(dag);
wolffd@0 100 msg = cell(1,N);
wolffd@0 101 observed = ~isemptycell(evidence(:));
wolffd@0 102
wolffd@0 103 for n=1:N
wolffd@0 104 ps = parents(dag, n);
wolffd@0 105 msg{n}.pi_from_parent = cell(1, length(ps));
wolffd@0 106 for i=1:length(ps)
wolffd@0 107 p = ps(i);
wolffd@0 108 msg{n}.pi_from_parent{i} = ones(ns(p), 1);
wolffd@0 109 end
wolffd@0 110
wolffd@0 111 cs = children(dag, n);
wolffd@0 112 msg{n}.lambda_from_child = cell(1, length(cs));
wolffd@0 113 for i=1:length(cs)
wolffd@0 114 c = cs(i);
wolffd@0 115 msg{n}.lambda_from_child{i} = ones(ns(n), 1);
wolffd@0 116 end
wolffd@0 117
wolffd@0 118 msg{n}.lambda = ones(ns(n), 1);
wolffd@0 119 msg{n}.pi = ones(ns(n), 1);
wolffd@0 120
wolffd@0 121 msg{n}.lambda_from_self = ones(ns(n), 1);
wolffd@0 122 end
wolffd@0 123
wolffd@0 124
wolffd@0 125 %%%%%%%%%
wolffd@0 126
wolffd@0 127 function msg = init_ev_msgs(engine, evidence, msg)
wolffd@0 128 % Initialize the lambdas with any evidence
wolffd@0 129
wolffd@0 130 [ss T] = size(evidence);
wolffd@0 131 bnet = bnet_from_engine(engine);
wolffd@0 132 pot_type = 'd';
wolffd@0 133 t = 1;
wolffd@0 134 hnodes = mysetdiff(1:ss, engine.onodes);
wolffd@0 135 for i=hnodes(:)'
wolffd@0 136 c = engine.obschild(i);
wolffd@0 137 if c > 0
wolffd@0 138 fam = family(bnet.dag, c);
wolffd@0 139 e = bnet.equiv_class(c, 1);
wolffd@0 140 CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,1));
wolffd@0 141 temp = pot_to_marginal(CPDpot);
wolffd@0 142 n = i;
wolffd@0 143 msg{n}.lambda_from_self = temp.T;
wolffd@0 144 end
wolffd@0 145 end
wolffd@0 146 for t=2:T
wolffd@0 147 for i=hnodes(:)'
wolffd@0 148 c = engine.obschild(i);
wolffd@0 149 if c > 0
wolffd@0 150 fam = family(bnet.dag, c, 2);
wolffd@0 151 e = bnet.equiv_class(c, 2);
wolffd@0 152 CPDpot = CPD_to_pot(pot_type, bnet.CPD{e}, fam, bnet.node_sizes(:), bnet.cnodes(:), evidence(:,t-1:t));
wolffd@0 153 temp = pot_to_marginal(CPDpot);
wolffd@0 154 n = i + (t-1)*ss;
wolffd@0 155 msg{n}.lambda_from_self = temp.T;
wolffd@0 156 end
wolffd@0 157 end
wolffd@0 158 end