nanstd.m

一个很有用的EM算法程序包

function sx = nanstd(x, flag)
%NANSTD   Standard deviation of available data, ignoring NaNs.
%
%    NANSTD(X) returns the standard deviation of the available data in
%    X, treating NaNs as missing values.  For vectors, NANSTD(X) is
%    the standard deviation of the non-NaN elements in X.  For
%    matrices, NANSTD(X) is a row vector containing the standard
%    deviation of the non-NaN elements in each column.
%
%    NANSTD(X) normalizes by (N-1) where, for each element of
%    NANSTD(X), N is number of available values.
%
%    NANSTD(X,0) normalizes by N and produces the second moment of the
%    available data about their mean.  NANSTD(X,1) is the same as
%    NANSTD(X).
%  
%    See also STD, NANMEAN.

  % maximum admissible fraction of missing values
  max_miss = 0.6;                

  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 nobs-1

  % if x is a vector, make sure it is a row vector
  if length(x)==prod(size(x))         
    x = x(:);                         
  end  
  [m,n]   = size(x);
  
  % replace NaNs with zeros.
  inan    = find(isnan(x));
  x(inan) = zeros(size(inan));
  
  % determine number of available observations on each variable
  [i,j]   = ind2sub([m,n], inan);     % subscripts of missing entries
  nans    = sparse(i,j,1,m,n);        % indicator matrix for missing values
  nobs    = m - sum(nans);

  % set nobs to NaN when there are too few entries to form robust average
  minobs  = m * (1 - max_miss);
  k       = find(nobs < minobs);
  nobs(k) = NaN;
  
  % center data
  xc      = x - repmat(sum(x) ./ nobs, m, 1);  % remove mean
  
  % standard deviation
  sx      = sqrt(sum(conj(xc).*xc) ./ (nobs-flag));