time   calls  line
                  1 function y = nanvar(x,w,dim)
                  2 %NANVAR Variance, ignoring NaNs.
                  3 %   Y = NANVAR(X) returns the sample variance of the values in X, treating
                  4 %   NaNs as missing values.  For a vector input, Y is the variance of the
                  5 %   non-NaN elements of X.  For a matrix input, Y is a row vector
                  6 %   containing the variance of the non-NaN elements in each column of X.
                  7 %   For N-D arrays, NANVAR operates along the first non-singleton dimension
                  8 %   of X.
                  9 %
                 10 %   NANVAR normalizes Y by N-1 if N>1, where N is the sample size of the 
                 11 %   non-NaN elements.  This is an unbiased estimator of the variance of the
                 12 %   population from which X is drawn, as long as X consists of independent,
                 13 %   identically distributed samples, and data are missing at random.  For
                 14 %   N=1, Y is normalized by N. 
                 15 %
                 16 %   Y = NANVAR(X,1) normalizes by N and produces the second moment of the
                 17 %   sample about its mean.  NANVAR(X,0) is the same as NANVAR(X).
                 18 %
                 19 %   Y = NANVAR(X,W) computes the variance using the weight vector W.  The
                 20 %   length of W must equal the length of the dimension over which NANVAR
                 21 %   operates, and its non-NaN elements must be nonnegative.  Elements of X
                 22 %   corresponding to NaN elements of W are ignored.
                 23 %
                 24 %   Y = NANVAR(X,W,DIM) takes the variance along dimension DIM of X.
                 25 %
                 26 %   See also VAR, NANSTD, NANMEAN, NANMEDIAN, NANMIN, NANMAX, NANSUM.
                 27 
                 28 %   Copyright 1984-2010 The MathWorks, Inc.
                 29 %   $Revision: 1.1.8.2 $  $Date: 2010/10/08 17:25:19 $
                 30 
             1   31 if nargin < 2 || isempty(w), w = 0; end 
                 32 
             1   33 sz = size(x); 
             1   34 if nargin < 3 || isempty(dim) 
                 35     % The output size for [] is a special case when DIM is not given.
             1   36     if isequal(x,[]), y = NaN(class(x)); return; end 
                 37 
                 38     % Figure out which dimension sum will work along.
             1   39     dim = find(sz ~= 1, 1); 
             1   40     if isempty(dim), dim = 1; end 
                 41 elseif dim > length(sz)
                 42     sz(end+1:dim) = 1;
                 43 end
                 44 
                 45 % Need to tile the mean of X to center it.
             1   46 tile = ones(size(sz)); 
             1   47 tile(dim) = sz(dim); 
                 48 
             1   49 if isequal(w,0) || isequal(w,1) 
                 50     % Count up non-NaNs.
             1   51     n = sum(~isnan(x),dim); 
                 52 
             1   53     if w == 0 
                 54         % The unbiased estimator: divide by (n-1).  Can't do this when
                 55         % n == 0 or 1, so n==1 => we'll return zeros
             1   56         denom = max(n-1, 1); 
                 57     else
                 58         % The biased estimator: divide by n.
                 59         denom = n; % n==1 => we'll return zeros
                 60     end
             1   61     denom(n==0) = NaN; % Make all NaNs return NaN, without a divideByZero warning 
                 62 
  0.01       1   63     x0 = x - repmat(nanmean(x, dim), tile); 
             1   64     y = nansum(abs(x0).^2, dim) ./ denom; % abs guarantees a real result 
                 65 
                 66 % Weighted variance
                 67 elseif numel(w) ~= sz(dim)
                 68     error(message('stats:nanvar:InvalidSizeWgts'));
                 69 elseif ~(isvector(w) && all(w(~isnan(w)) >= 0))
                 70     error(message('stats:nanvar:InvalidWgts'));
                 71 else
                 72     % Embed W in the right number of dims.  Then replicate it out along the
                 73     % non-working dims to match X's size.
                 74     wresize = ones(size(sz)); wresize(dim) = sz(dim);
                 75     wtile = sz; wtile(dim) = 1;
                 76     w = repmat(reshape(w, wresize), wtile);
                 77 
                 78     % Count up non-NaNs.
                 79     n = nansum(~isnan(x).*w,dim);
                 80 
                 81     x0 = x - repmat(nansum(w.*x, dim) ./ n, tile);
                 82     y = nansum(w .* abs(x0).^2, dim) ./ n; % abs guarantees a real result
                 83 end