|
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 Gaussian node
|
|
wolffd@0
|
3 % function CPD = update_ess(CPD, fmarginal, evidence, ns, cnodes, hidden_bitv)
|
|
wolffd@0
|
4
|
|
wolffd@0
|
5 %if nargin < 6
|
|
wolffd@0
|
6 % hidden_bitv = zeros(1, max(fmarginal.domain));
|
|
wolffd@0
|
7 % hidden_bitv(find(isempty(evidence)))=1;
|
|
wolffd@0
|
8 %end
|
|
wolffd@0
|
9
|
|
wolffd@0
|
10 dom = fmarginal.domain;
|
|
wolffd@0
|
11 self = dom(end);
|
|
wolffd@0
|
12 ps = dom(1:end-1);
|
|
wolffd@0
|
13 cps = myintersect(ps, cnodes);
|
|
wolffd@0
|
14 dps = mysetdiff(ps, cps);
|
|
wolffd@0
|
15
|
|
wolffd@0
|
16 CPD.nsamples = CPD.nsamples + 1;
|
|
wolffd@0
|
17 [ss cpsz dpsz] = size(CPD.weights); % ss = self size
|
|
wolffd@0
|
18 [ss dpsz] = size(CPD.mean);
|
|
wolffd@0
|
19
|
|
wolffd@0
|
20 % Let X be the cts parent (if any), Y be the cts child (self).
|
|
wolffd@0
|
21
|
|
wolffd@0
|
22 if ~hidden_bitv(self) & ~any(hidden_bitv(cps)) & all(hidden_bitv(dps))
|
|
wolffd@0
|
23 % Speedup for the common case that all cts nodes are observed, all discrete nodes are hidden
|
|
wolffd@0
|
24 % Since X and Y are observed, SYY = 0, SXX = 0, SXY = 0
|
|
wolffd@0
|
25 % Since discrete parents are hidden, we do not need to add evidence to w.
|
|
wolffd@0
|
26 w = fmarginal.T(:);
|
|
wolffd@0
|
27 CPD.Wsum = CPD.Wsum + w;
|
|
wolffd@0
|
28 y = evidence{self};
|
|
wolffd@0
|
29 Cyy = y*y';
|
|
wolffd@0
|
30 if ~CPD.useC
|
|
wolffd@0
|
31 WY = repmat(w(:)',ss,1); % WY(y,i) = w(i)
|
|
wolffd@0
|
32 WYY = repmat(reshape(WY, [ss 1 dpsz]), [1 ss 1]); % WYY(y,y',i) = w(i)
|
|
wolffd@0
|
33 %CPD.WYsum = CPD.WYsum + WY .* repmat(y(:), 1, dpsz);
|
|
wolffd@0
|
34 CPD.WYsum = CPD.WYsum + y(:) * w(:)';
|
|
wolffd@0
|
35 CPD.WYYsum = CPD.WYYsum + WYY .* repmat(reshape(Cyy, [ss ss 1]), [1 1 dpsz]);
|
|
wolffd@0
|
36 else
|
|
wolffd@0
|
37 W = w(:)';
|
|
wolffd@0
|
38 W2 = reshape(W, [1 1 dpsz]);
|
|
wolffd@0
|
39 CPD.WYsum = CPD.WYsum + rep_mult(W, y(:), size(CPD.WYsum));
|
|
wolffd@0
|
40 CPD.WYYsum = CPD.WYYsum + rep_mult(W2, Cyy, size(CPD.WYYsum));
|
|
wolffd@0
|
41 end
|
|
wolffd@0
|
42 if cpsz > 0 % X exists
|
|
wolffd@0
|
43 x = cat(1, evidence{cps}); x = x(:);
|
|
wolffd@0
|
44 Cxx = x*x';
|
|
wolffd@0
|
45 Cxy = x*y';
|
|
wolffd@0
|
46 WX = repmat(w(:)',cpsz,1); % WX(x,i) = w(i)
|
|
wolffd@0
|
47 WXX = repmat(reshape(WX, [cpsz 1 dpsz]), [1 cpsz 1]); % WXX(x,x',i) = w(i)
|
|
wolffd@0
|
48 WXY = repmat(reshape(WX, [cpsz 1 dpsz]), [1 ss 1]); % WXY(x,y,i) = w(i)
|
|
wolffd@0
|
49 if ~CPD.useC
|
|
wolffd@0
|
50 CPD.WXsum = CPD.WXsum + WX .* repmat(x(:), 1, dpsz);
|
|
wolffd@0
|
51 CPD.WXXsum = CPD.WXXsum + WXX .* repmat(reshape(Cxx, [cpsz cpsz 1]), [1 1 dpsz]);
|
|
wolffd@0
|
52 CPD.WXYsum = CPD.WXYsum + WXY .* repmat(reshape(Cxy, [cpsz ss 1]), [1 1 dpsz]);
|
|
wolffd@0
|
53 else
|
|
wolffd@0
|
54 CPD.WXsum = CPD.WXsum + rep_mult(W, x(:), size(CPD.WXsum));
|
|
wolffd@0
|
55 CPD.WXXsum = CPD.WXXsum + rep_mult(W2, Cxx, size(CPD.WXXsum));
|
|
wolffd@0
|
56 CPD.WXYsum = CPD.WXYsum + rep_mult(W2, Cxy, size(CPD.WXYsum));
|
|
wolffd@0
|
57 end
|
|
wolffd@0
|
58 end
|
|
wolffd@0
|
59 return;
|
|
wolffd@0
|
60 end
|
|
wolffd@0
|
61
|
|
wolffd@0
|
62 % general (non-vectorized) case
|
|
wolffd@0
|
63 fullm = add_evidence_to_gmarginal(fmarginal, evidence, ns, cnodes); % slow!
|
|
wolffd@0
|
64
|
|
wolffd@0
|
65 if dpsz == 1 % no discrete parents
|
|
wolffd@0
|
66 w = 1;
|
|
wolffd@0
|
67 else
|
|
wolffd@0
|
68 w = fullm.T(:);
|
|
wolffd@0
|
69 end
|
|
wolffd@0
|
70
|
|
wolffd@0
|
71 CPD.Wsum = CPD.Wsum + w;
|
|
wolffd@0
|
72 xi = 1:cpsz;
|
|
wolffd@0
|
73 yi = (cpsz+1):(cpsz+ss);
|
|
wolffd@0
|
74 for i=1:dpsz
|
|
wolffd@0
|
75 muY = fullm.mu(yi, i);
|
|
wolffd@0
|
76 SYY = fullm.Sigma(yi, yi, i);
|
|
wolffd@0
|
77 CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY;
|
|
wolffd@0
|
78 CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y]
|
|
wolffd@0
|
79 if cpsz > 0
|
|
wolffd@0
|
80 muX = fullm.mu(xi, i);
|
|
wolffd@0
|
81 SXX = fullm.Sigma(xi, xi, i);
|
|
wolffd@0
|
82 SXY = fullm.Sigma(xi, yi, i);
|
|
wolffd@0
|
83 CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX;
|
|
wolffd@0
|
84 CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX');
|
|
wolffd@0
|
85 CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY');
|
|
wolffd@0
|
86 end
|
|
wolffd@0
|
87 end
|
|
wolffd@0
|
88
|
