wolffd@0: function [engine, loglik] = enter_evidence(engine, evidence, varargin)
wolffd@0: % ENTER_EVIDENCE Add the specified evidence to the network (gaussian_inf_engine)
wolffd@0: % [engine, loglik] = enter_evidence(engine, evidence, ...)
wolffd@0: %
wolffd@0: % evidence{i} = [] if if X(i) is hidden, and otherwise contains its observed value (scalar or column vector)
wolffd@0: 
wolffd@0: bnet = bnet_from_engine(engine);
wolffd@0: ns = bnet.node_sizes;
wolffd@0: O = find(~isemptycell(evidence));
wolffd@0: H = find(isemptycell(evidence));
wolffd@0: vals = cat(1, evidence{O});
wolffd@0: 
wolffd@0: % Compute Pr(H|o)
wolffd@0: [Hmu, HSigma, loglik] = condition_gaussian(engine.mu, engine.Sigma, H, O, vals(:), ns);
wolffd@0: 
wolffd@0: engine.Hmu = Hmu;
wolffd@0: engine.HSigma = HSigma;
wolffd@0: engine.hnodes = H;
wolffd@0: 
wolffd@0: %%%%%%%%
wolffd@0: 
wolffd@0: function [mu2, Sigma2, loglik] = condition_gaussian(mu, Sigma, X, Y, y, ns)
wolffd@0: % CONDITION_GAUSSIAN Compute Pr(X|Y=y) where X and Y are jointly Gaussian.
wolffd@0: % [mu2, Sigma2, ll] = condition_gaussian(mu, Sigma, X, Y, y, ns)
wolffd@0: 
wolffd@0: if isempty(y)
wolffd@0:   mu2 = mu;
wolffd@0:   Sigma2 = Sigma;
wolffd@0:   loglik = 0;
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: use_log = 1;
wolffd@0: 
wolffd@0: if length(Y)==length(mu) % instantiating every variable
wolffd@0:   mu2 = y;
wolffd@0:   Sigma2 = zeros(length(y));
wolffd@0:   loglik = gaussian_prob(y, mu, Sigma, use_log);
wolffd@0:   return;
wolffd@0: end
wolffd@0: 
wolffd@0: [muX, muY, SXX, SXY, SYX, SYY] = partition_matrix_vec(mu, Sigma, X, Y, ns);
wolffd@0: K = SXY*inv(SYY);
wolffd@0: mu2 = muX + K*(y-muY);
wolffd@0: Sigma2 = SXX - K*SYX;
wolffd@0: loglik = gaussian_prob(y, muY, SYY, use_log);