Covariance reconstruction

function demoReconCov
%reconstruct covariance matrix with unobserved elements
    rng(42) %seed
    mu      = [1 2 3 4];
    NUN_smp = 10000;        % sample size
    for c=1:3 %loop through the cases
        %% Define generating true model
        matC_True{1} = [1.5 -0.3 -0.23 0.2;
                        -0.3 0.99 0.11 -0.04;
                        -0.23 0.11 1.3 -0.14;
                        0.2 -0.04 -0.14 1.09;];
        %% Reconstruction using Wishart priors by sega_sai @https://stats.stackexchange.com/q/371743
        matC_Wish{1} = [1.501 -0.299 0.126 -0.311;
                        -0.299 0.988 0.112 -0.039;
                        0.126 0.112 1.301 -0.138;
                        -0.311 -0.039 -0.138 1.09];

        matC_True{2} = [1.19 -0.3 -0.23 0.2;
                        -0.3 0.68 0.11 -0.04;
                        -0.23 0.11 0.9 -0.14;
                        0.2 -0.04 -0.14 0.65];

        matC_Wish{2} = [1.189 -0.305 -0.201 -0.294;
                        -0.305 0.684 0.107 -0.046;
                        -0.201 0.107 0.901 -0.141;
                        -0.294 -0.046 -0.141 0.653];

        matC_True{3} = [1.19 -0.3 -0.23 0.2;
                        -0.3 0.68 0.11 -0.04;
                        -0.23 0.11 0.9 -0.14;
                        0.2 -0.04 -0.14 0.65];

        matC_Wish{3} = [1.198 -0.312 0.599 0.402;
                        -0.312 0.682 0.119 -0.035;
                        0.599 0.119 0.9 -0.145;
                        0.402 -0.035 -0.145 0.647];

        %% Experiment #1
        x12     = mvnrnd(mu,matC_True{c},NUN_smp);
        x12     = x12(:,1:2);  % only x_1, x_2 are observed
        mu12    = mean(x12);   % estimate mvn model
        matC12  = cov(x12);
        %% Experiment #2
        x234    = mvnrnd(mu,matC_True{c},NUN_smp);
        x234    = x234(:,2:4); % only x_2, x_3, x_4 are observed
        mu234   = mean(x234);  % estimate mvn model
        matC234 = cov(x234);
        %% Determinant maximization
        cellCovMu   = {matC12,mu12;...  % cell array with Covariance and Mu of each experiment
                       matC234,mu234};
        cellID      = {[1 2]';[2 3 4]'}; % observed dimensions in each experiment
        [corMatNew,muVec] = maxDetMvn(cellCovMu,cellID); %determinant maximization

        %% Validation, all x observed, calculate likelihoods
        xVal        = mvnrnd(mu,matC_True{c},100);
        logLTrue    = sum(log(mvnpdf(xVal,mu,matC_True{c})));           % true model likelihood
        logLEstWish = sum(log(mvnpdf(xVal,muVec',matC_Wish{c})));    % Method #1: wishart prior
        logLEstMaxD = sum(log(mvnpdf(xVal,muVec',corMatNew)));  % Method #2: max determinant
        disp(['=== Case ' num2str(c) ' ====']);
        disp(['True model logL:      ' num2str(logLTrue,'%.2f')] );
        disp(['Wishart prior logL:   ' num2str(logLEstWish,'%.2f') ]);
        disp(['max Determinant logL: ' num2str(logLEstMaxD,'%.2f') ]);
        disp(['']);
    end
%%
    function [nlogD, corMat] = logDetCor(corMat,vIJ,x)
    %calculates negative log determinant of a correlation matrix
    %substituting elements i=vIJ(:,1),j=vIJ(:,1) with x(:)
        tmpIdx  = sub2ind(size(corMat),...
                    [vIJ(:,1); vIJ(:,2)],...
                    [vIJ(:,2); vIJ(:,1)]);
        corMat(tmpIdx)=[x;x];
        eigV    = eig(corMat);
        if(any(eigV<0))
            nlogD = 10000;
        else
            nlogD = -log(prod(eigV));
        end
    end
%%
    function [corMatNew,muVec] = maxDetMvn(cellCovMu,cellID)
    %determinant maximization
    %cellCovMu: cell array containing covariance and mean vector of each
    %           measurement
    %cellID:    cell array containing obveserved dimensions in each measurement

        NUM_cv      = size(cellCovMu,1);
        MAX_pr      = max(cell2mat(cellID));
        cellElmCov  = cell(NUM_cv,2);
        cellElmMu   = cell(NUM_cv,2);
        for i=1:NUM_cv
            [ti,tj] = find(~isnan(cellCovMu{i,1}));
            tmpPrID = cellID{i};
            tmpIdx  = sub2ind([1 1]*MAX_pr,tmpPrID(ti),tmpPrID(tj));
            cellElmCov{i,1} = tmpIdx;
            cellElmCov{i,2} = cellCovMu{i,1}(:);
            cellElmMu{i,1}  = tmpPrID(:);
            cellElmMu{i,2}  = cellCovMu{i,2}(:);
        end

        covMat  = nan(MAX_pr,MAX_pr);
        vecCidx = cell2mat(cellElmCov(:,1));
        vecCval = cell2mat(cellElmCov(:,2));
        covMat  = meanElm(vecCidx,vecCval,covMat); %average overlapping elements

        muVec   = nan(MAX_pr,1);
        vecMidx = cell2mat(cellElmMu(:,1));
        vecMval = cell2mat(cellElmMu(:,2));
        muVec   = meanElm(vecMidx,vecMval,muVec);

        [vSig,corMat]   = cov2corr(covMat);
        [vi,vj]         = find(isnan(corMat));
        keepI   = vi>vj;
        numX    = sum(keepI);   %number of unknowns
        vIJ     = [vi(keepI) vj(keepI)];
        lb  = -1*ones(numX,1);  %lower bound
        ub  = ones(numX,1);     %upper bound
        x0  = ones(numX,1);     %initial guess
        fun = @(x)logDetCor(corMat,vIJ,x);
        opts= optimoptions('fmincon','Display','off');
        x   = fmincon(fun,x0,[],[],[],[],lb,ub,[],opts);   %constrained minimization
        [~,corMatNew] = logDetCor(corMat,vIJ,x);        %max entropy correlation
        corMatNew     = corr2cov(vSig,corMatNew);       %convert back to covariance
        end
%%
    function inMat = meanElm(vecIdx,vecVal,inMat)
    %average repeating matrix elements
        [uI,~,aI] = unique(vecIdx);
        numU = histc(vecIdx,[uI-0.5; uI(end)+0.5]);
        for i=1:numel(uI)
            if(numU(i)>1)
                inMat(uI(i)) = sum(vecVal(aI==i))/numU(i);
            else
                inMat(uI(i)) = vecVal(aI==i);
            end
        end
    end
end