wolffd@0: % Test gmux.
wolffd@0: % The following model, where Y is a gmux node,
wolffd@0: % and M is set to 1, should be equivalent to X1 -> Y
wolffd@0: %
wolffd@0: % X1 Xn M
wolffd@0: % \ |  /
wolffd@0: %   Y
wolffd@0: 
wolffd@0: n = 3;
wolffd@0: N = n+2;
wolffd@0: Xs = 1:n;
wolffd@0: M = n+1; 
wolffd@0: Y = n+2;
wolffd@0: dag = zeros(N,N);
wolffd@0: dag([Xs M], Y)=1; 
wolffd@0: 
wolffd@0: dnodes = M;
wolffd@0: ns = zeros(1, N);
wolffd@0: sz = 2;
wolffd@0: ns(Xs) = sz;
wolffd@0: ns(M) = n;
wolffd@0: ns(Y) = sz;
wolffd@0: 
wolffd@0: bnet = mk_bnet(dag, ns, 'discrete', M, 'observed', [M Y]);
wolffd@0: 
wolffd@0: psz = ns(Xs(1));
wolffd@0: selfsz = ns(Y);
wolffd@0: 
wolffd@0: W = randn(selfsz, psz);
wolffd@0: mu = randn(selfsz, 1);
wolffd@0: Sigma = eye(selfsz, selfsz);
wolffd@0: 
wolffd@0: bnet.CPD{M} = root_CPD(bnet, M);
wolffd@0: for i=Xs(:)'
wolffd@0:   bnet.CPD{i} = gaussian_CPD(bnet, i, 'mean', zeros(psz, 1), 'cov', eye(psz, psz));
wolffd@0: end
wolffd@0: bnet.CPD{Y} = gmux_CPD(bnet, Y, 'mean', mu, 'weights', W, 'cov', Sigma);
wolffd@0:   
wolffd@0: evidence = cell(1,N);
wolffd@0: yval = randn(selfsz, 1);
wolffd@0: evidence{Y} = yval;
wolffd@0: m = 2;
wolffd@0: %notm = not(m-1)+1; % only valid for n=2
wolffd@0: notm = mysetdiff(1:n, m);
wolffd@0: evidence{M} = m;
wolffd@0: 
wolffd@0: engines = {};
wolffd@0: engines{end+1} = jtree_inf_engine(bnet);
wolffd@0: engines{end+1} = pearl_inf_engine(bnet, 'protocol', 'parallel');
wolffd@0: 
wolffd@0: for e=1:length(engines)
wolffd@0:   engines{e} = enter_evidence(engines{e}, evidence);
wolffd@0:   mXm{e} = marginal_nodes(engines{e}, Xs(m));
wolffd@0: 
wolffd@0:   % Since M=m, only Xm was updated.
wolffd@0:   % Hence the posterior on Xnotm should equal the prior.
wolffd@0:   for i=notm(:)'
wolffd@0:     mXnotm = marginal_nodes(engines{e}, Xs(i));
wolffd@0:     assert(approxeq(mXnotm.mu, zeros(psz,1)))
wolffd@0:     assert(approxeq(mXnotm.Sigma, eye(psz, psz)))
wolffd@0:   end
wolffd@0: end
wolffd@0: 
wolffd@0: % Check that all engines give the same posterior
wolffd@0: for e=2:length(engines)
wolffd@0:   assert(approxeq(mXm{e}.mu, mXm{1}.mu))
wolffd@0:   assert(approxeq(mXm{e}.Sigma, mXm{1}.Sigma))
wolffd@0: end
wolffd@0: 
wolffd@0: 
wolffd@0: % Compute the correct posterior by building Xm -> Y
wolffd@0: 
wolffd@0: N = 2;
wolffd@0: dag = zeros(N,N);
wolffd@0: dag(1, 2)=1;
wolffd@0: ns = [psz selfsz];
wolffd@0: bnet = mk_bnet(dag, ns, 'discrete', [], 'observed', 2);
wolffd@0: 
wolffd@0: bnet.CPD{1} = gaussian_CPD(bnet, 1, 'mean', zeros(psz, 1), 'cov', eye(psz, psz));
wolffd@0: bnet.CPD{2} = gaussian_CPD(bnet, 2, 'mean', mu, 'cov', Sigma, 'weights', W);
wolffd@0: 
wolffd@0: jengine  = jtree_inf_engine(bnet);
wolffd@0: evidence = {[], yval};
wolffd@0: jengine = enter_evidence(jengine, evidence); % apply Bayes rule to invert the arc
wolffd@0: mX = marginal_nodes(jengine, 1);
wolffd@0: 
wolffd@0: for e=1:length(engines)
wolffd@0:   assert(approxeq(mX.mu, mXm{e}.mu))
wolffd@0:   assert(approxeq(mX.Sigma, mXm{e}.Sigma))
wolffd@0: end