Mercurial > hg > camir-aes2014

annotate toolboxes/FullBNT-1.0.7/bnt/CPDs/@mlp_CPD/update_ess.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 CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
wolffd@0 2 % UPDATE_ESS Update the Expected Sufficient Statistics of a CPD (MLP)
wolffd@0 3 % CPD = update_ess(CPD, family_marginal, evidence, node_sizes, cnodes, hidden_bitv)
wolffd@0 4 %
wolffd@0 5 % fmarginal = overall posterior distribution of self and its parents
wolffd@0 6 % fmarginal(i1,i2...,ik,s)=prob(Pa1=i1,...,Pak=ik, self=s| X)
wolffd@0 7 %
wolffd@0 8 % => 1) prob(self|Pa1,...,Pak)=fmarginal/prob(Pa1,...,Pak) with prob(Pa1,...,Pak)=sum{s,fmarginal}
wolffd@0 9 % [self estimation -> CPD.self_vals]
wolffd@0 10 % 2) prob(Pa1,...,Pak) [SCG weights -> CPD.eso_weights]
wolffd@0 11 %
wolffd@0 12 % Hidden_bitv is ignored
wolffd@0 13
wolffd@0 14 % Written by Pierpaolo Brutti
wolffd@0 15
wolffd@0 16 if ~adjustable_CPD(CPD), return; end
wolffd@0 17
wolffd@0 18 dom = fmarginal.domain;
wolffd@0 19 cdom = myintersect(dom, cnodes);
wolffd@0 20 assert(~any(isemptycell(evidence(cdom))));
wolffd@0 21 ns(cdom)=1;
wolffd@0 22
wolffd@0 23 self = dom(end);
wolffd@0 24 ps=dom(1:end-1);
wolffd@0 25 dpdom=mysetdiff(ps,cdom);
wolffd@0 26
wolffd@0 27 dnodes = mysetdiff(1:length(ns), cnodes);
wolffd@0 28
wolffd@0 29 ddom = myintersect(ps, dnodes); %
wolffd@0 30 if isempty(evidence{self}), % if self is hidden in what follow we must
wolffd@0 31 ddom = myintersect(dom, dnodes); % consider its dimension
wolffd@0 32 end %
wolffd@0 33
wolffd@0 34 odom = dom(~isemptycell(evidence(dom)));
wolffd@0 35 hdom = dom(isemptycell(evidence(dom))); % hidden parents in domain
wolffd@0 36
wolffd@0 37 dobs = myintersect(ddom, odom);
wolffd@0 38 dvals = cat(1, evidence{dobs});
wolffd@0 39 ens = ns; % effective node sizes
wolffd@0 40 ens(dobs) = 1;
wolffd@0 41
wolffd@0 42 dpsz=prod(ns(dpdom));
wolffd@0 43 S=prod(ens(ddom));
wolffd@0 44 subs = ind2subv(ens(ddom), 1:S);
wolffd@0 45 mask = find_equiv_posns(dobs, ddom);
wolffd@0 46 for i=1:length(mask),
wolffd@0 47 subs(:,mask(i)) = dvals(i);
wolffd@0 48 end
wolffd@0 49 supportedQs = subv2ind(ns(ddom), subs);
wolffd@0 50
wolffd@0 51 Qarity = prod(ns(ddom));
wolffd@0 52 if isempty(ddom),
wolffd@0 53 Qarity = 1;
wolffd@0 54 end
wolffd@0 55 fullm.T = zeros(Qarity, 1);
wolffd@0 56 fullm.T(supportedQs) = fmarginal.T(:);
wolffd@0 57
wolffd@0 58 % For dynamic (recurrent) net-------------------------------------------------------------
wolffd@0 59 % ----------------------------------------------------------------------------------------
wolffd@0 60 high=size(evidence,1); % slice height
wolffd@0 61 ss_ns=ns(1:high); % single slice nodes sizes
wolffd@0 62 pos=self; %
wolffd@0 63 slice_num=0; %
wolffd@0 64 while pos>high, %
wolffd@0 65 slice_num=slice_num+1; % find active slice
wolffd@0 66 pos=pos-high; % pos=self posistion into a single slice
wolffd@0 67 end %
wolffd@0 68
wolffd@0 69 last_dim=pos-1; %
wolffd@0 70 if isempty(evidence{self}), %
wolffd@0 71 last_dim=pos; %
wolffd@0 72 end % last_dim=last reshaping dimension
wolffd@0 73 reg=dom-slice_num*high;
wolffd@0 74 dex=myintersect(reg(find(reg>=0)), [1:last_dim]); %
wolffd@0 75 rs_dim=ss_ns(dex); % reshaping dimensions
wolffd@0 76
wolffd@0 77 if slice_num>0,
wolffd@0 78 act_slice=[]; past_ancest=[]; %
wolffd@0 79 act_slice=slice_num*high+[1:high]; % recover the active slice nodes
wolffd@0 80 % past_ancest=mysetdiff(ddom, act_slice);
wolffd@0 81 past_ancest=mysetdiff(ps, act_slice); % recover ancestors contained into past slices
wolffd@0 82 app=ns(past_ancest);
wolffd@0 83 rs_dim=[app(:)' rs_dim(:)']; %
wolffd@0 84 end %
wolffd@0 85 if length(rs_dim)==1, rs_dim=[1 rs_dim]; end %
wolffd@0 86 if size(rs_dim,1)~=1, rs_dim=rs_dim'; end %
wolffd@0 87
wolffd@0 88 fullm.T=reshape(fullm.T, rs_dim); % reshaping the marginal
wolffd@0 89
wolffd@0 90 % ----------------------------------------------------------------------------------------
wolffd@0 91 % ----------------------------------------------------------------------------------------
wolffd@0 92
wolffd@0 93 % X = cts parent, R = discrete self
wolffd@0 94
wolffd@0 95 % 1) observations vector -> CPD.parents_vals -------------------------------------------------
wolffd@0 96 x = cat(1, evidence{cdom});
wolffd@0 97
wolffd@0 98 % 2) weights vector -> CPD.eso_weights -------------------------------------------------------
wolffd@0 99 if isempty(evidence{self}) % R is hidden
wolffd@0 100 sum_over=length(rs_dim);
wolffd@0 101 app=sum(fullm.T, sum_over);
wolffd@0 102 pesi=reshape(app,[dpsz,1]);
wolffd@0 103 clear app;
wolffd@0 104 else
wolffd@0 105 pesi=reshape(fullm.T,[dpsz,1]);
wolffd@0 106 end
wolffd@0 107
wolffd@0 108 assert(approxeq(sum(pesi),1));
wolffd@0 109
wolffd@0 110 % 3) estimate (if R is hidden) or recover (if R is obs) self'value----------------------------
wolffd@0 111 if isempty(evidence{self}) % R is hidden
wolffd@0 112 app=mk_stochastic(fullm.T); % P(self|Pa1,...,Pak)=fmarginal/prob(Pa1,...,Pak)
wolffd@0 113 app=reshape(app,[dpsz ns(self)]); % matrix size: prod{j,ns(Paj)} x ns(self)
wolffd@0 114 r=app;
wolffd@0 115 clear app;
wolffd@0 116 else
wolffd@0 117 r = zeros(dpsz,ns(self));
wolffd@0 118 for i=1:dpsz
wolffd@0 119 if pesi(i)~=0, r(i,evidence{self}) = 1; end
wolffd@0 120 end
wolffd@0 121 end
wolffd@0 122 for i=1:dpsz
wolffd@0 123 if pesi(i) ~=0, assert(approxeq(sum(r(i,:)),1)); end
wolffd@0 124 end
wolffd@0 125
wolffd@0 126 CPD.nsamples = CPD.nsamples + 1;
wolffd@0 127 CPD.parent_vals(CPD.nsamples,:) = x(:)';
wolffd@0 128 for i=1:dpsz
wolffd@0 129 CPD.eso_weights(CPD.nsamples,:,i)=pesi(i);
wolffd@0 130 CPD.self_vals(CPD.nsamples,:,i) = r(i,:);
wolffd@0 131 end