function [mu, Sigma, ll] = naive_gaussian_marginal(bnet, nodes)
% NAIVE_GAUSSIAN_MARGINAL Compute the global multivariate Gaussian corresponding to bnet and marginalize it.
% [mu, Sigma, ll] = naive_gaussian_marginal(bnet, nodes)
% 
% ll is the log-likelihood of the data.
% This is the Gaussian version of naive_discrete_marginal (but much faster!).

num_nodes = length(bnet.adj_mat);
all_nodes = 1:num_nodes;
assert(subset(bnet.cts_nodes, all_nodes));
ns = bnet.node_sizes;
obs = bnet.values(bnet.observed_nodes);

H = bnet.hidden_nodes;
O = bnet.observed_nodes;
% Compute Pr(H,O)
mvg = gbn_to_mvg(bnet);
% Compute Pr(H|o)
[Hmu, HSigma, ll] = condition_gaussian(mvg.mu, mvg.Sigma, H, O, all_nodes, obs, ns);
% Compute Pr(Hkeep|o)
Hkeep = intersect(H, nodes);
Hsum = setdiff(H, Hkeep);
[ll2, mu, Sigma] = marginalize_gaussian(ll, Hmu, HSigma, Hkeep, Hsum, H, ns(H));
% Add Okeep to Pr(Hkeep|o)
Okeep = intersect(O, nodes);
Okeep_vals = bnet.values(Okeep);
[mu, Sigma] = add_ev_to_gaussian(mu, Sigma, Hkeep, Okeep, nodes, Okeep_vals, ns(nodes));

