|
wolffd@0
|
1 % Test gmux.
|
|
wolffd@0
|
2 % The following model, where Y is a gmux node,
|
|
wolffd@0
|
3 % and M is set to 1, should be equivalent to X1 -> Y
|
|
wolffd@0
|
4 %
|
|
wolffd@0
|
5 % X1 Xn M
|
|
wolffd@0
|
6 % \ | /
|
|
wolffd@0
|
7 % Y
|
|
wolffd@0
|
8
|
|
wolffd@0
|
9 n = 3;
|
|
wolffd@0
|
10 N = n+2;
|
|
wolffd@0
|
11 Xs = 1:n;
|
|
wolffd@0
|
12 M = n+1;
|
|
wolffd@0
|
13 Y = n+2;
|
|
wolffd@0
|
14 dag = zeros(N,N);
|
|
wolffd@0
|
15 dag([Xs M], Y)=1;
|
|
wolffd@0
|
16
|
|
wolffd@0
|
17 dnodes = M;
|
|
wolffd@0
|
18 ns = zeros(1, N);
|
|
wolffd@0
|
19 sz = 2;
|
|
wolffd@0
|
20 ns(Xs) = sz;
|
|
wolffd@0
|
21 ns(M) = n;
|
|
wolffd@0
|
22 ns(Y) = sz;
|
|
wolffd@0
|
23
|
|
wolffd@0
|
24 bnet = mk_bnet(dag, ns, 'discrete', M, 'observed', [M Y]);
|
|
wolffd@0
|
25
|
|
wolffd@0
|
26 psz = ns(Xs(1));
|
|
wolffd@0
|
27 selfsz = ns(Y);
|
|
wolffd@0
|
28
|
|
wolffd@0
|
29 W = randn(selfsz, psz);
|
|
wolffd@0
|
30 mu = randn(selfsz, 1);
|
|
wolffd@0
|
31 Sigma = eye(selfsz, selfsz);
|
|
wolffd@0
|
32
|
|
wolffd@0
|
33 bnet.CPD{M} = root_CPD(bnet, M);
|
|
wolffd@0
|
34 for i=Xs(:)'
|
|
wolffd@0
|
35 bnet.CPD{i} = gaussian_CPD(bnet, i, 'mean', zeros(psz, 1), 'cov', eye(psz, psz));
|
|
wolffd@0
|
36 end
|
|
wolffd@0
|
37 bnet.CPD{Y} = gmux_CPD(bnet, Y, 'mean', mu, 'weights', W, 'cov', Sigma);
|
|
wolffd@0
|
38
|
|
wolffd@0
|
39 evidence = cell(1,N);
|
|
wolffd@0
|
40 yval = randn(selfsz, 1);
|
|
wolffd@0
|
41 evidence{Y} = yval;
|
|
wolffd@0
|
42 m = 2;
|
|
wolffd@0
|
43 %notm = not(m-1)+1; % only valid for n=2
|
|
wolffd@0
|
44 notm = mysetdiff(1:n, m);
|
|
wolffd@0
|
45 evidence{M} = m;
|
|
wolffd@0
|
46
|
|
wolffd@0
|
47 engines = {};
|
|
wolffd@0
|
48 engines{end+1} = jtree_inf_engine(bnet);
|
|
wolffd@0
|
49 engines{end+1} = pearl_inf_engine(bnet, 'protocol', 'parallel');
|
|
wolffd@0
|
50
|
|
wolffd@0
|
51 for e=1:length(engines)
|
|
wolffd@0
|
52 engines{e} = enter_evidence(engines{e}, evidence);
|
|
wolffd@0
|
53 mXm{e} = marginal_nodes(engines{e}, Xs(m));
|
|
wolffd@0
|
54
|
|
wolffd@0
|
55 % Since M=m, only Xm was updated.
|
|
wolffd@0
|
56 % Hence the posterior on Xnotm should equal the prior.
|
|
wolffd@0
|
57 for i=notm(:)'
|
|
wolffd@0
|
58 mXnotm = marginal_nodes(engines{e}, Xs(i));
|
|
wolffd@0
|
59 assert(approxeq(mXnotm.mu, zeros(psz,1)))
|
|
wolffd@0
|
60 assert(approxeq(mXnotm.Sigma, eye(psz, psz)))
|
|
wolffd@0
|
61 end
|
|
wolffd@0
|
62 end
|
|
wolffd@0
|
63
|
|
wolffd@0
|
64 % Check that all engines give the same posterior
|
|
wolffd@0
|
65 for e=2:length(engines)
|
|
wolffd@0
|
66 assert(approxeq(mXm{e}.mu, mXm{1}.mu))
|
|
wolffd@0
|
67 assert(approxeq(mXm{e}.Sigma, mXm{1}.Sigma))
|
|
wolffd@0
|
68 end
|
|
wolffd@0
|
69
|
|
wolffd@0
|
70
|
|
wolffd@0
|
71 % Compute the correct posterior by building Xm -> Y
|
|
wolffd@0
|
72
|
|
wolffd@0
|
73 N = 2;
|
|
wolffd@0
|
74 dag = zeros(N,N);
|
|
wolffd@0
|
75 dag(1, 2)=1;
|
|
wolffd@0
|
76 ns = [psz selfsz];
|
|
wolffd@0
|
77 bnet = mk_bnet(dag, ns, 'discrete', [], 'observed', 2);
|
|
wolffd@0
|
78
|
|
wolffd@0
|
79 bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', zeros(psz, 1), 'cov', eye(psz, psz));
|
|
wolffd@0
|
80 bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', mu, 'cov', Sigma, 'weights', W);
|
|
wolffd@0
|
81
|
|
wolffd@0
|
82 jengine = jtree_inf_engine(bnet);
|
|
wolffd@0
|
83 evidence = {[], yval};
|
|
wolffd@0
|
84 jengine = enter_evidence(jengine, evidence); % apply Bayes rule to invert the arc
|
|
wolffd@0
|
85 mX = marginal_nodes(jengine, 1);
|
|
wolffd@0
|
86
|
|
wolffd@0
|
87 for e=1:length(engines)
|
|
wolffd@0
|
88 assert(approxeq(mX.mu, mXm{e}.mu))
|
|
wolffd@0
|
89 assert(approxeq(mX.Sigma, mXm{e}.Sigma))
|
|
wolffd@0
|
90 end
|
