Joint MVN

%% Define generating true model
rng(1) %seed
mu      = [1 2 3 4];
matC    = [1.0 0.2 0.1 0.2;
           0.2 0.5 0.3 0.2;
           0.1 0.3 0.3 0.1;
           0.2 0.2 0.1 1.0;];
%% Experiement #1
x12     = mvnrnd(mu,matC,10000);
x12     = x12(:,1:2);  % only x_1, x_2 are observed
mu12    = mean(x12);   % estimate mvn model
matC12  = cov(x12);
%% Experiement #2
x234    = mvnrnd(mu,matC,10000);
x234    = x234(:,2:4); % only x_2, x_3, x_4 are observed
mu234   = mean(x234);  % estimate mvn model
matC234 = cov(x234);

%% Estimate joint probability
P1  = [1 0 0 0; 0 1 0 0];           %projection matrices
P2  = [0 1 0 0; 0 0 1 0; 0 0 0 1];

matCEst     = inv(...
                  ((P1')*inv(matC12)*P1) + ...
                  ((P2')*inv(matC234)*P2)...
                  );
muEst       = matCEst * ...
                  ((P1')*inv(matC12)*mu12' + ...
                   (P2')*inv(matC234)*mu234');

%% Validation, all x observed, calculate likelihoods
xVal        = mvnrnd(mu,matC,100);
logLTrue    = sum(log(mvnpdf(xVal,mu,matC)));        % -417.6038
logLEst     = sum(log(mvnpdf(xVal,muEst',matCEst))); % -435.6095