|
wolffd@0
|
1 function smallpot = marginalize_pot(bigpot, keep)
|
|
wolffd@0
|
2 % MARGINALIZE_POT Marginalize a cgpot onto a smaller domain.
|
|
wolffd@0
|
3 % smallpot = marginalize_pot(bigpot, keep)
|
|
wolffd@0
|
4
|
|
wolffd@0
|
5 sumover = mysetdiff(bigpot.domain, keep);
|
|
wolffd@0
|
6 csumover = myintersect(sumover, bigpot.cdom);
|
|
wolffd@0
|
7 dsumover = myintersect(sumover, bigpot.ddom);
|
|
wolffd@0
|
8 dkeep = myintersect(keep, bigpot.ddom);
|
|
wolffd@0
|
9 ckeep = myintersect(keep, bigpot.cdom);
|
|
wolffd@0
|
10 %ns = sparse(1, max(bigpot.domain)); % must be full, so I is an integer
|
|
wolffd@0
|
11 ns = zeros(1, max(bigpot.domain));
|
|
wolffd@0
|
12 ns(bigpot.ddom) = bigpot.dsizes;
|
|
wolffd@0
|
13 ns(bigpot.cdom) = bigpot.csizes;
|
|
wolffd@0
|
14
|
|
wolffd@0
|
15 % sum(ns(csumover))==0 is like isempty(csumover) but handles observed nodes.
|
|
wolffd@0
|
16 % Similarly, prod(ns(dsumover))==1 is like isempty(dsumover)
|
|
wolffd@0
|
17
|
|
wolffd@0
|
18 % Marginalize the cts parts.
|
|
wolffd@0
|
19 % If we are in canonical form, we stay that way, since moment form might not exist.
|
|
wolffd@0
|
20 % Besides, we would like to minimize the number of conversions.
|
|
wolffd@0
|
21 if sum(ns(csumover)) > 0
|
|
wolffd@0
|
22 if bigpot.subtype == 'm'
|
|
wolffd@0
|
23 for i=1:bigpot.dsize
|
|
wolffd@0
|
24 bigpot.mom{i} = marginalize_pot(bigpot.mom{i}, ckeep);
|
|
wolffd@0
|
25 end
|
|
wolffd@0
|
26 else
|
|
wolffd@0
|
27 for i=1:bigpot.dsize
|
|
wolffd@0
|
28 bigpot.can{i} = marginalize_pot(bigpot.can{i}, ckeep);
|
|
wolffd@0
|
29 end
|
|
wolffd@0
|
30 end
|
|
wolffd@0
|
31 end
|
|
wolffd@0
|
32
|
|
wolffd@0
|
33 % If we are not marginalizing over any discrete nodes, we are done.
|
|
wolffd@0
|
34 if prod(ns(dsumover))==1
|
|
wolffd@0
|
35 smallpot = cgpot(dkeep, ckeep, ns, bigpot.can, bigpot.mom, bigpot.subtype);
|
|
wolffd@0
|
36 return;
|
|
wolffd@0
|
37 end
|
|
wolffd@0
|
38
|
|
wolffd@0
|
39 % To marginalize the discrete parts, we must be in moment form.
|
|
wolffd@0
|
40 bigpot = cg_can_to_mom(bigpot);
|
|
wolffd@0
|
41
|
|
wolffd@0
|
42 I = prod(ns(dkeep));
|
|
wolffd@0
|
43 J = prod(ns(dsumover));
|
|
wolffd@0
|
44 C = sum(ns(ckeep));
|
|
wolffd@0
|
45
|
|
wolffd@0
|
46 % Reshape bigpot into the form mu1(:,j,i), where i is in dkeep, j is in dsumover
|
|
wolffd@0
|
47 T1 = zeros(I,J);
|
|
wolffd@0
|
48 mu1 = zeros(C,J,I);
|
|
wolffd@0
|
49 Sigma1 = zeros(C,C,J,I);
|
|
wolffd@0
|
50 sum_map = find_equiv_posns(dsumover, bigpot.ddom);
|
|
wolffd@0
|
51 keep_map = find_equiv_posns(dkeep, bigpot.ddom);
|
|
wolffd@0
|
52 iv = zeros(1, length(bigpot.ddom)); % index vector
|
|
wolffd@0
|
53 for i=1:I
|
|
wolffd@0
|
54 keep_iv = ind2subv(ns(dkeep), i);
|
|
wolffd@0
|
55 iv(keep_map) = keep_iv;
|
|
wolffd@0
|
56 for j=1:J
|
|
wolffd@0
|
57 sum_iv = ind2subv(ns(dsumover), j);
|
|
wolffd@0
|
58 iv(sum_map) = sum_iv;
|
|
wolffd@0
|
59 k = subv2ind(ns(bigpot.ddom), iv);
|
|
wolffd@0
|
60 mom = struct(bigpot.mom{k}); % violate object privacy
|
|
wolffd@0
|
61 T1(i,j) = exp(mom.logp);
|
|
wolffd@0
|
62 if C > 0 % so mu1 and Sigma1 are non-empty
|
|
wolffd@0
|
63 mu1(:,j,i) = mom.mu;
|
|
wolffd@0
|
64 Sigma1(:,:,j,i) = mom.Sigma;
|
|
wolffd@0
|
65 end
|
|
wolffd@0
|
66 end
|
|
wolffd@0
|
67 end
|
|
wolffd@0
|
68
|
|
wolffd@0
|
69 % Collapse the mixture of Gaussians
|
|
wolffd@0
|
70 coef = mk_stochastic(T1); % coef must be convex combination
|
|
wolffd@0
|
71 T2 = sum(T1,2);
|
|
wolffd@0
|
72 T2 = T2 + (T2==0)*eps;
|
|
wolffd@0
|
73 %if C > 0, disp('collapsing onto '); disp(leep); end
|
|
wolffd@0
|
74 mu = [];
|
|
wolffd@0
|
75 Sigma = [];
|
|
wolffd@0
|
76 mom = cell(1,I);
|
|
wolffd@0
|
77 for i=1:I
|
|
wolffd@0
|
78 if C > 0
|
|
wolffd@0
|
79 [mu, Sigma] = collapse_mog(mu1(:,:,i), Sigma1(:,:,:,i), coef(i,:));
|
|
wolffd@0
|
80 end
|
|
wolffd@0
|
81 logp = log(T2(i));
|
|
wolffd@0
|
82 mom{i} = mpot(ckeep, ns(ckeep), logp, mu, Sigma);
|
|
wolffd@0
|
83 end
|
|
wolffd@0
|
84
|
|
wolffd@0
|
85 smallpot = cgpot(dkeep, ckeep, ns, [], mom, 'm');
|
|
wolffd@0
|
86
|
