function [c, xm, xc] = nancov(x, flag);%NANCOV  Sample covariance matrix of available data, ignoring NaNs.% %    NANCOV(X) returns the covariance matrix of X estimated from all%    available data in X, ignoring NaN-elements that indicate missing%    values. If X is a vector, NANCOV(X) returns the variance. For%    matrices, where each row is an observation and each column a%    variable, NANCOV(X) is the covariance matrix estimate from all%    available data.%%    NANCOV(X) normalizes each covariance matrix element by (Nt-1)%    where Nt is the number of available product terms contributing to%    the covariance matrix element.%%    NANCOV(X,0) normalizes by Nt and produces the second moment%    matrix of the observations about their mean.  NANCOV(X,1) is the%    same as NANCOV(X).%  %    If all terms in the sum of the product and/or cross-product terms%    contributing to a covariance matrix element are missing, the%    covariance matrix element is set to zero.%%    Because the covariance matrix elements are normalized with%    different normalization factors, depending on how many terms%    contribute to them, the covariance matrix is not necessarily%    positive semidefinite.%%    [C, M, XC] = NANCOV(X) returns the covariance matrix estimate%    C from all available data, the mean M of the available data,%    and the centered data matrix XC with the mean removed from%    each column. %%    See also COV, CENTER, NANSTD, NANMEAN, NANSUM.     error(nargchk(1,2,nargin))          % check number of input arguments   if isempty(x)                       % check for empty input.    sx = NaN;    return  end  if ndims(x) > 2,  error('Data must be vector or 2-D array.'); end    if nargin < 2, flag = 1; end        % default: normalize by nt-1  if length(x)==prod(size(x))         % if x is vector, make sure    x = x(:);                         % it is a row vector  end     % center data  [xc,xm] = center(x);  [m,n]   = size(xc);    % replace NaNs with zeros.  inan    = find(isnan(xc));  xc(inan)= zeros(size(inan));      nonnan  = ones(m,n);                % indicator matrix showing where there  nonnan(inan) = zeros(size(inan));   % are non-missing  elements in xc  nt      = nonnan' * nonnan - flag;  % normalization factor for each 				      % covariance matrix element      % set covariance to 0 when there are no terms to form sum  k       = find(nt < 1);  nt(k)   = 1;  c       = xc' * xc ./ nt;