function [mu, K_tot] = multivarSampleStats(x, K_x, nosparse)
%x is a DxN set of N measurements of a D dimensional random vector
%K_x is a cell array of the inverse covariances of each measurement
%mu is the weighted sample mean of these measurements
%K_mu is the weighted sample inverse covariance of these measurements

if ~exist('nosparse', 'var')
    nosparse = false;
end

if issparse(K_x{1})
    for cK = 1:length(K_x)
        K_x{cK} = full(K_x{cK});
    end
end

good_measurements = true(1,length(K_x));
for ii = 1:length(K_x)
    if any(isnan(K_x{ii}(:))) || any(isnan(x(:,ii)))
        good_measurements(ii) = false;
    end
end

if ~any(good_measurements)
    mu = nan(size(x,1), 1);
    K_tot = nan(size(x,1));
    return;
end



K_x = K_x(good_measurements);
x = x(:,good_measurements);



K_x = cell2mat(reshape(K_x,1,1,[]));

K_tot = sum(K_x,3);


mu = zeros(size(x,1),1);
for ii = 1:size(x,2)
    mu = mu + K_x(:,:,ii)*x(:,ii);
end

mu = K_tot\mu;

C = zeros(size(mu,1));
for ii = 1:size(x,2) 
    C = C + K_x(:,:,ii)*(x(:,ii) - mu)*(x(:,ii) - mu)';
end
C = inv(K_tot) + K_tot\C;
K_mu = inv(C);

if ~nosparse
    K_mu = sparse(K_mu);
elseif nosparse
    K_mu = full(K_mu);
end

end