function [muXY, SigmaXY] = condGaussToJoint(muX, SigmaX, muY, SigmaY, Wgivc)% Compute P(X,Y) from P(X) * P(Y|X) where P(X)=N(X;muX,SigmaX) % and P(Y|X) = N(Y; WX + muY, SigmaY)% For details on how to compute a Gaussian from a Bayes net% - "Gaussian Influence Diagrams", R. Shachter and C. R. Kenley, Management Science, 35(5):527--550, 1989.% size(W) = dy x dxdx = length(muX);dy = length(muY);muXY = [muX(:); muY];WXY = [zeros(dx,dx) W';     zeros(dx,dx) zeros(dy,dy)];D = [SigmaX zeros(dy,dy);	   zeros(dx,dx) SigmaY];U = inv(eye(size(W)) - W')';Sigma = U' * D * U;