octave_packages/statistics-1.1.3/nanstd.m

   1 ## Copyright (C) 2001 Paul Kienzle <pkienzle@users.sf.net>
   2 ##
   3 ## This program is free software; you can redistribute it and/or modify it under
   4 ## the terms of the GNU General Public License as published by the Free Software
   5 ## Foundation; either version 3 of the License, or (at your option) any later
   6 ## version.
   7 ##
   8 ## This program is distributed in the hope that it will be useful, but WITHOUT
   9 ## ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
  10 ## FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public License for more
  11 ## details.
  12 ##
  13 ## You should have received a copy of the GNU General Public License along with
  14 ## this program; if not, see <http://www.gnu.org/licenses/>.
  15
  16 ## -*- texinfo -*-
  17 ## @deftypefn {Function File} {@var{v} =} nanstd (@var{X})
  18 ## @deftypefnx{Function File} {@var{v} =} nanstd (@var{X}, @var{opt})
  19 ## @deftypefnx{Function File} {@var{v} =} nanstd (@var{X}, @var{opt}, @var{dim})
  20 ## Compute the standard deviation while ignoring NaN values.
  21 ##
  22 ## @code{nanstd} is identical to the @code{std} function except that NaN values are
  23 ## ignored.  If all values are NaN, the standard deviation is returned as NaN.
  24 ## If there is only a single non-NaN value, the deviation is returned as 0.
  25 ##
  26 ## The argument @var{opt} determines the type of normalization to use. Valid values
  27 ## are
  28 ##
  29 ## @table @asis
  30 ## @item 0:
  31 ##   normalizes with @math{N-1}, provides the square root of best unbiased estimator of
  32 ##   the variance [default]
  33 ## @item 1:
  34 ##   normalizes with @math{N}, this provides the square root of the second moment around
  35 ##   the mean
  36 ## @end table
  37 ##
  38 ## The third argument @var{dim} determines the dimension along which the standard
  39 ## deviation is calculated.
  40 ##
  41 ## @seealso{std, nanmin, nanmax, nansum, nanmedian, nanmean}
  42 ## @end deftypefn
  43
  44 function v = nanstd (X, opt, varargin)
  45   if nargin < 1
  46     print_usage;
  47   else
  48     if nargin < 3
  49       dim = min(find(size(X)>1));
  50       if isempty(dim), dim=1; endif;
  51     else
  52       dim = varargin{1};
  53     endif
  54     if ((nargin < 2) || isempty(opt))
  55       opt = 0;
  56     endif
  57
  58     ## determine the number of non-missing points in each data set
  59     n = sum (!isnan(X), varargin{:});
  60
  61     ## replace missing data with zero and compute the mean
  62     X(isnan(X)) = 0;
  63     meanX = sum (X, varargin{:}) ./ n;
  64
  65     ## subtract the mean from the data and compute the sum squared
  66     sz = ones(1,length(size(X)));
  67     sz(dim) = size(X,dim);
  68     v = sumsq (X - repmat(meanX,sz), varargin{:});
  69
  70     ## because the missing data was set to zero each missing data
  71     ## point will contribute (-meanX)^2 to sumsq, so remove these
  72     v = v - (meanX .^ 2) .* (size(X,dim) - n);
  73
  74     if (opt == 0)
  75       ## compute the standard deviation from the corrected sumsq using
  76       ## max(n-1,1) in the denominator so that the std for a single point is 0
  77       v = sqrt ( v ./ max(n - 1, 1) );
  78     elseif (opt == 1)
  79       ## compute the standard deviation from the corrected sumsq
  80       v = sqrt ( v ./ n );
  81     else
  82       error ("std: unrecognized normalization type");
  83     endif
  84
  85   endif
  86 endfunction