function C = covnan(A) % same as cov, but leaving out NaNs. % divides by zero, if matrix has too many missing values. [s, n] = sumnan(A, 1); B = A - ones(size(A, 1), 1) * (s ./ n); nb = size(B, 2); C = zeros(nb); for q = 1 : nb, [s, n] = sumnan(B(:, q : end) .* repmat(B(:, q), 1, nb - q + 1), 1); C(q, q : end) = s ./ n; % was (n-1) C(q : end, q) = C(q, q : end)'; % use symmetry for optimization end;