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