1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
|
function CPD = update_tied_ess(CPD, domain, engine, evidence, ns, cnodes)
if ~adjustable_CPD(CPD), return; end
nCPDs = size(domain, 2);
fmarginal = cell(1, nCPDs);
for l=1:nCPDs
fmarginal{l} = marginal_family(engine, nodes(l));
end
[ss cpsz dpsz] = size(CPD.weights);
if const_evidence_pattern(engine)
dom = domain(:,1);
dnodes = mysetdiff(1:length(ns), cnodes);
ddom = myintersect(dom, dnodes);
cdom = myintersect(dom, cnodes);
odom = dom(~isemptycell(evidence(dom)));
hdom = dom(isemptycell(evidence(dom)));
% If all hidden nodes are discrete and all cts nodes are observed
% (e.g., HMM with Gaussian output)
% we can add the observed evidence in parallel
if mysubset(ddom, hdom) & mysubset(cdom, odom)
[mu, Sigma, T] = add_cts_ev_to_marginals(fmarginal, evidence, ns, cnodes);
else
mu = zeros(ss, dpsz, nCPDs);
Sigma = zeros(ss, ss, dpsz, nCPDs);
T = zeros(dpsz, nCPDs);
for l=1:nCPDs
[mu(:,:,l), Sigma(:,:,:,l), T(:,l)] = add_ev_to_marginals(fmarginal{l}, evidence, ns, cnodes);
end
end
end
CPD.nsamples = CPD.nsamples + nCPDs;
if dpsz == 1 % no discrete parents
w = 1;
else
w = fullm.T(:);
end
CPD.Wsum = CPD.Wsum + w;
% Let X be the cts parent (if any), Y be the cts child (self).
xi = 1:cpsz;
yi = (cpsz+1):(cpsz+ss);
for i=1:dpsz
muY = fullm.mu(yi, i);
SYY = fullm.Sigma(yi, yi, i);
CPD.WYsum(:,i) = CPD.WYsum(:,i) + w(i)*muY;
CPD.WYYsum(:,:,i) = CPD.WYYsum(:,:,i) + w(i)*(SYY + muY*muY'); % E[X Y] = Cov[X,Y] + E[X] E[Y]
if cpsz > 0
muX = fullm.mu(xi, i);
SXX = fullm.Sigma(xi, xi, i);
SXY = fullm.Sigma(xi, yi, i);
CPD.WXsum(:,i) = CPD.WXsum(:,i) + w(i)*muX;
CPD.WXYsum(:,:,i) = CPD.WXYsum(:,:,i) + w(i)*(SXY + muX*muY');
CPD.WXXsum(:,:,i) = CPD.WXXsum(:,:,i) + w(i)*(SXX + muX*muX');
end
end
%%%%%%%%%%%%%
function fullm = add_evidence_to_marginal(fmarginal, evidence, ns, cnodes)
dom = fmarginal.domain;
% Find out which values of the discrete parents (if any) are compatible with
% the discrete evidence (if any).
dnodes = mysetdiff(1:length(ns), cnodes);
ddom = myintersect(dom, dnodes);
cdom = myintersect(dom, cnodes);
odom = dom(~isemptycell(evidence(dom)));
hdom = dom(isemptycell(evidence(dom)));
dobs = myintersect(ddom, odom);
dvals = cat(1, evidence{dobs});
ens = ns; % effective node sizes
ens(dobs) = 1;
S = prod(ens(ddom));
subs = ind2subv(ens(ddom), 1:S);
mask = find_equiv_posns(dobs, ddom);
subs(mask) = dvals;
supportedQs = subv2ind(ns(ddom), subs);
if isempty(ddom)
Qarity = 1;
else
Qarity = prod(ns(ddom));
end
fullm.T = zeros(Qarity, 1);
fullm.T(supportedQs) = fmarginal.T(:);
% Now put the hidden cts parts into their right blocks,
% leaving the observed cts parts as 0.
cobs = myintersect(cdom, odom);
chid = myintersect(cdom, hdom);
cvals = cat(1, evidence{cobs});
n = sum(ns(cdom));
fullm.mu = zeros(n,Qarity);
fullm.Sigma = zeros(n,n,Qarity);
if ~isempty(chid)
chid_blocks = block(find_equiv_posns(chid, cdom), ns(cdom));
end
if ~isempty(cobs)
cobs_blocks = block(find_equiv_posns(cobs, cdom), ns(cdom));
end
for i=1:length(supportedQs)
Q = supportedQs(i);
if ~isempty(chid)
fullm.mu(chid_blocks, Q) = fmarginal.mu(:, i);
fullm.Sigma(chid_blocks, chid_blocks, Q) = fmarginal.Sigma(:,:,i);
end
if ~isempty(cobs)
fullm.mu(cobs_blocks, Q) = cvals(:);
end
end
|