covnan.m

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;