SHARE
TWEET

Covariance reconstruction

a guest Oct 19th, 2018 11 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. function demoReconCov
  2. %reconstruct covariance matrix with unobserved elements
  3.     rng(42) %seed
  4.     mu      = [1 2 3 4];
  5.     NUN_smp = 10000;        % sample size
  6.     for c=1:3 %loop through the cases
  7.         %% Define generating true model
  8.         matC_True{1} = [1.5 -0.3 -0.23 0.2;
  9.                         -0.3 0.99 0.11 -0.04;
  10.                         -0.23 0.11 1.3 -0.14;
  11.                         0.2 -0.04 -0.14 1.09;];
  12.         %% Reconstruction using Wishart priors by sega_sai @https://stats.stackexchange.com/q/371743
  13.         matC_Wish{1} = [1.501 -0.299 0.126 -0.311;
  14.                         -0.299 0.988 0.112 -0.039;
  15.                         0.126 0.112 1.301 -0.138;
  16.                         -0.311 -0.039 -0.138 1.09];
  17.  
  18.         matC_True{2} = [1.19 -0.3 -0.23 0.2;
  19.                         -0.3 0.68 0.11 -0.04;
  20.                         -0.23 0.11 0.9 -0.14;
  21.                         0.2 -0.04 -0.14 0.65];
  22.  
  23.         matC_Wish{2} = [1.189 -0.305 -0.201 -0.294;
  24.                         -0.305 0.684 0.107 -0.046;
  25.                         -0.201 0.107 0.901 -0.141;
  26.                         -0.294 -0.046 -0.141 0.653];
  27.  
  28.         matC_True{3} = [1.19 -0.3 -0.23 0.2;
  29.                         -0.3 0.68 0.11 -0.04;
  30.                         -0.23 0.11 0.9 -0.14;
  31.                         0.2 -0.04 -0.14 0.65];
  32.  
  33.         matC_Wish{3} = [1.198 -0.312 0.599 0.402;
  34.                         -0.312 0.682 0.119 -0.035;
  35.                         0.599 0.119 0.9 -0.145;
  36.                         0.402 -0.035 -0.145 0.647];
  37.        
  38.         %% Experiment #1
  39.         x12     = mvnrnd(mu,matC_True{c},NUN_smp);
  40.         x12     = x12(:,1:2);  % only x_1, x_2 are observed
  41.         mu12    = mean(x12);   % estimate mvn model
  42.         matC12  = cov(x12);
  43.         %% Experiment #2
  44.         x234    = mvnrnd(mu,matC_True{c},NUN_smp);
  45.         x234    = x234(:,2:4); % only x_2, x_3, x_4 are observed
  46.         mu234   = mean(x234);  % estimate mvn model
  47.         matC234 = cov(x234);
  48.         %% Determinant maximization
  49.         cellCovMu   = {matC12,mu12;...  % cell array with Covariance and Mu of each experiment
  50.                        matC234,mu234};
  51.         cellID      = {[1 2]';[2 3 4]'}; % observed dimensions in each experiment
  52.         [corMatNew,muVec] = maxDetMvn(cellCovMu,cellID); %determinant maximization
  53.  
  54.         %% Validation, all x observed, calculate likelihoods
  55.         xVal        = mvnrnd(mu,matC_True{c},100);
  56.         logLTrue    = sum(log(mvnpdf(xVal,mu,matC_True{c})));           % true model likelihood
  57.         logLEstWish = sum(log(mvnpdf(xVal,muVec',matC_Wish{c})));    % Method #1: wishart prior
  58.         logLEstMaxD = sum(log(mvnpdf(xVal,muVec',corMatNew)));  % Method #2: max determinant
  59.         disp(['=== Case ' num2str(c) ' ====']);
  60.         disp(['True model logL:      ' num2str(logLTrue,'%.2f')] );
  61.         disp(['Wishart prior logL:   ' num2str(logLEstWish,'%.2f') ]);
  62.         disp(['max Determinant logL: ' num2str(logLEstMaxD,'%.2f') ]);
  63.         disp(['']);
  64.     end
  65. %%
  66.     function [nlogD, corMat] = logDetCor(corMat,vIJ,x)
  67.     %calculates negative log determinant of a correlation matrix
  68.     %substituting elements i=vIJ(:,1),j=vIJ(:,1) with x(:)
  69.         tmpIdx  = sub2ind(size(corMat),...
  70.                     [vIJ(:,1); vIJ(:,2)],...
  71.                     [vIJ(:,2); vIJ(:,1)]);
  72.         corMat(tmpIdx)=[x;x];
  73.         eigV    = eig(corMat);
  74.         if(any(eigV<0))
  75.             nlogD = 10000;
  76.         else
  77.             nlogD = -log(prod(eigV));
  78.         end
  79.     end
  80. %%
  81.     function [corMatNew,muVec] = maxDetMvn(cellCovMu,cellID)
  82.     %determinant maximization
  83.     %cellCovMu: cell array containing covariance and mean vector of each
  84.     %           measurement
  85.     %cellID:    cell array containing obveserved dimensions in each measurement
  86.  
  87.         NUM_cv      = size(cellCovMu,1);
  88.         MAX_pr      = max(cell2mat(cellID));
  89.         cellElmCov  = cell(NUM_cv,2);
  90.         cellElmMu   = cell(NUM_cv,2);
  91.         for i=1:NUM_cv
  92.             [ti,tj] = find(~isnan(cellCovMu{i,1}));
  93.             tmpPrID = cellID{i};
  94.             tmpIdx  = sub2ind([1 1]*MAX_pr,tmpPrID(ti),tmpPrID(tj));
  95.             cellElmCov{i,1} = tmpIdx;
  96.             cellElmCov{i,2} = cellCovMu{i,1}(:);
  97.             cellElmMu{i,1}  = tmpPrID(:);
  98.             cellElmMu{i,2}  = cellCovMu{i,2}(:);
  99.         end
  100.  
  101.         covMat  = nan(MAX_pr,MAX_pr);
  102.         vecCidx = cell2mat(cellElmCov(:,1));
  103.         vecCval = cell2mat(cellElmCov(:,2));
  104.         covMat  = meanElm(vecCidx,vecCval,covMat); %average overlapping elements
  105.  
  106.         muVec   = nan(MAX_pr,1);
  107.         vecMidx = cell2mat(cellElmMu(:,1));
  108.         vecMval = cell2mat(cellElmMu(:,2));
  109.         muVec   = meanElm(vecMidx,vecMval,muVec);
  110.  
  111.         [vSig,corMat]   = cov2corr(covMat);    
  112.         [vi,vj]         = find(isnan(corMat));
  113.         keepI   = vi>vj;
  114.         numX    = sum(keepI);   %number of unknowns
  115.         vIJ     = [vi(keepI) vj(keepI)];
  116.         lb  = -1*ones(numX,1);  %lower bound
  117.         ub  = ones(numX,1);     %upper bound
  118.         x0  = ones(numX,1);     %initial guess
  119.         fun = @(x)logDetCor(corMat,vIJ,x);
  120.         opts= optimoptions('fmincon','Display','off');
  121.         x   = fmincon(fun,x0,[],[],[],[],lb,ub,[],opts);   %constrained minimization
  122.         [~,corMatNew] = logDetCor(corMat,vIJ,x);        %max entropy correlation
  123.         corMatNew     = corr2cov(vSig,corMatNew);       %convert back to covariance
  124.         end
  125. %%
  126.     function inMat = meanElm(vecIdx,vecVal,inMat)
  127.     %average repeating matrix elements
  128.         [uI,~,aI] = unique(vecIdx);
  129.         numU = histc(vecIdx,[uI-0.5; uI(end)+0.5]);
  130.         for i=1:numel(uI)
  131.             if(numU(i)>1)
  132.                 inMat(uI(i)) = sum(vecVal(aI==i))/numU(i);
  133.             else
  134.                 inMat(uI(i)) = vecVal(aI==i);
  135.             end
  136.         end
  137.     end
  138. end
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top