octave_packages/nan-2.5.5/moment.m

   1 function M=moment(i,p,opt,DIM)
   2 % MOMENT estimates the p-th moment
   3 %
   4 % M = moment(x, p [,opt] [,DIM])
   5 % M = moment(H, p [,opt])
   6 %   calculates p-th central moment from data x in dimension DIM
   7 %       of from Histogram H
   8 %
   9 % p     moment of order p
  10 % opt   'ac': absolute 'a' and/or central ('c') moment
  11 %       DEFAULT: '' raw moments are estimated
  12 % DIM   dimension
  13 %       1: STATS of columns
  14 %       2: STATS of rows
  15 %       default or []: first DIMENSION, with more than 1 element
  16 %
  17 % features:
  18 % - can deal with NaN's (missing values)
  19 % - dimension argument
  20 % - compatible to Matlab and Octave
  21 %
  22 % see also: STD, VAR, SKEWNESS, KURTOSIS, STATISTIC,
  23 %
  24 % REFERENCE(S):
  25 % http://mathworld.wolfram.com/Moment.html
  26
  27 %    $Id: moment.m 8223 2011-04-20 09:16:06Z schloegl $
  28 %    Copyright (C) 2000-2002,2010 by Alois Schloegl <alois.schloegl@gmail.com>
  29 %    This functions is part of the NaN-toolbox
  30 %    http://pub.ist.ac.at/~schloegl/matlab/NaN/
  31
  32 %    This program is free software; you can redistribute it and/or modify
  33 %    it under the terms of the GNU General Public License as published by
  34 %    the Free Software Foundation; either version 3 of the License, or
  35 %    (at your option) any later version.
  36 %
  37 %    This program is distributed in the hope that it will be useful,
  38 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  39 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  40 %    GNU General Public License for more details.
  41 %
  42 %    You should have received a copy of the GNU General Public License
  43 %    along with this program; If not, see <http://www.gnu.org/licenses/>.
  44
  45 if nargin==2,
  46         DIM=[];
  47         opt=[];
  48 elseif nargin==3,
  49         DIM=[];
  50 elseif nargin==4,
  51
  52 else
  53         fprintf('Error MOMENT: invalid number of arguments\n');
  54         return;
  55 end;
  56
  57 if p<=0;
  58         fprintf('Error MOMENT: invalid model order p=%f\n',p);
  59         return;
  60 end;
  61
  62 if isnumeric(opt) || ~isnumeric(DIM),
  63         tmp = DIM;
  64         DIM = opt;
  65         opt = tmp;
  66 end;
  67 if isempty(opt),
  68         opt='r';
  69 end;
  70 if isempty(DIM),
  71         DIM = find(size(i)>1,1);
  72         if isempty(DIM), DIM=1; end;
  73 end;
  74
  75 N = nan;
  76 if isstruct(i),
  77     if isfield(i,'HISTOGRAM'),
  78         sz = size(i.H)./size(i.X);
  79         X  = repmat(i.X,sz);
  80         if any(opt=='c'),
  81                 N = sumskipnan(i.H,1);  % N
  82                 N = max(N-1,0);         % for unbiased estimation
  83                 S = sumskipnan(i.H.*X,1);       % sum
  84                 X = X - repmat(S./N, size(X)./size(S)); % remove mean
  85         end;
  86         if any(opt=='a'),
  87                 X = abs(X);
  88         end;
  89         [M,n] = sumskipnan(X.^p.*i.H,1);
  90     else
  91         warning('invalid datatype')
  92     end;
  93 else
  94         if any(opt=='c'),
  95                 [S,N] = sumskipnan(i,DIM);      % gemerate N and SUM
  96                 N = max(N-1,0);                 % for unbiased estimation
  97                 i = i - repmat(S./N, size(i)./size(S)); % remove mean
  98         end;
  99         if any(opt=='a'),
 100                 i = abs(i);
 101         end;
 102         [M,n] = sumskipnan(i.^p,DIM);
 103 end;
 104
 105 if isnan(N), N=n; end;
 106 M = M./N;