octave_packages/nan-2.5.5/statistic.m

   1 function [varargout]=statistic(i,DIM,fun)
   2 % STATISTIC estimates various statistics at once.
   3 %
   4 % R = STATISTIC(x,DIM)
   5 %   calculates all statistic (see list of fun) in dimension DIM
   6 %   R is a struct with all statistics
   7 %
   8 % y = STATISTIC(x,fun)
   9 %   estimate of fun on dimension DIM
  10 %   y gives the statistic of fun
  11 %
  12 % DIM   dimension
  13 %       1: STATS of columns
  14 %       2: STATS of rows
  15 %       N: STATS of  N-th dimension
  16 %       default or []: first DIMENSION, with more than 1 element
  17 %
  18 % fun   'mean'  mean
  19 %       'std'   standard deviation
  20 %       'var'   variance
  21 %       'sem'   standard error of the mean
  22 %       'rms'   root mean square
  23 %       'meansq' mean of squares
  24 %       'sum'   sum
  25 %       'sumsq' sum of squares
  26 %       'CM#'   central moment of order #
  27 %       'skewness' skewness
  28 %       'kurtosis' excess coefficient (Fisher kurtosis)
  29 %       'mad'   mean absolute deviation
  30 %
  31 % features:
  32 % - can deal with NaN's (missing values)
  33 % - dimension argument
  34 % - compatible to Matlab and Octave
  35 %
  36 % see also: SUMSKIPNAN
  37 %
  38 % REFERENCE(S):
  39 % [1] http://www.itl.nist.gov/
  40 % [2] http://mathworld.wolfram.com/
  41
  42 %       $Id: statistic.m 8223 2011-04-20 09:16:06Z schloegl $
  43 %       Copyright (C) 2000-2003,2010 by Alois Schloegl <alois.schloegl@gmail.com>
  44 %       This function is part of the NaN-toolbox
  45 %       http://pub.ist.ac.at/~schloegl/matlab/NaN/
  46
  47 %    This program is free software; you can redistribute it and/or modify
  48 %    it under the terms of the GNU General Public License as published by
  49 %    the Free Software Foundation; either version 3 of the License, or
  50 %    (at your option) any later version.
  51 %
  52 %    This program is distributed in the hope that it will be useful,
  53 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  54 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  55 %    GNU General Public License for more details.
  56 %
  57 %    You should have received a copy of the GNU General Public License
  58 %    along with this program; If not, see <http://www.gnu.org/licenses/>.
  59
  60
  61
  62 if nargin==1,
  63         DIM=[];
  64         fun=[];
  65 elseif nargin==2,
  66         if ~isnumeric(DIM),
  67                 fun=DIM;
  68                 DIM=[];
  69         else
  70                 fun=[];
  71         end
  72 end
  73 if isempty(DIM),
  74         DIM = find(size(i)>1,1);
  75         if isempty(DIM), DIM=1; end;
  76 end;
  77
  78 %R.N    = sumskipnan(~isnan(i),DIM);    % number of elements
  79 [R.SUM,R.N,R.SSQ] = sumskipnan(i,DIM);  % sum
  80 %R.S3P  = sumskipnan(i.^3,DIM);         % sum of 3rd power
  81 R.S4P  = sumskipnan(i.^4,DIM);          % sum of 4th power
  82 %R.S5P  = sumskipnan(i.^5,DIM);         % sum of 5th power
  83
  84 R.MEAN  = R.SUM./R.N;                   % mean
  85 R.MSQ   = R.SSQ./R.N;                   % mean square
  86 R.RMS   = sqrt(R.MSQ);                  % root mean square
  87 %R.SSQ0 = R.SSQ-R.SUM.*R.MEAN;          % sum square of mean removed
  88 R.SSQ0  = R.SSQ - real(R.SUM).*real(R.MEAN) - imag(R.SUM).*imag(R.MEAN);        % sum square of mean removed
  89
  90 %if flag_implicit_unbiased_estim;    %% ------- unbiased estimates -----------
  91     n1  = max(R.N-1,0);                 % in case of n=0 and n=1, the (biased) variance, STD and SEM are INF
  92 %else
  93 %    n1 = R.N;
  94 %end;
  95
  96 R.VAR   = R.SSQ0./n1;                   % variance (unbiased)
  97 R.STD   = sqrt(R.VAR);                  % standard deviation
  98 R.SEM   = sqrt(R.SSQ0./(R.N.*n1));      % standard error of the mean
  99 R.SEV   = sqrt(n1.*(n1.*R.S4P./R.N+(R.N.^2-2*R.N+3).*(R.SSQ./R.N).^2)./(R.N.^3)); % standard error of the variance
 100 R.COEFFICIENT_OF_VARIATION = R.STD./R.MEAN;
 101
 102 q = quantile(i, (1:3)/4, DIM);
 103
 104 %sz=size(i);sz(DIM)=1;
 105 %Q0500=repmat(nan,sz);
 106 %Q0250=Q0500;
 107 %Q0750=Q0500;
 108 %MODE=Q0500;
 109 %for k=1:size(i,2),
 110 %        tmp = sort(i(:,k));
 111         %ix = find(~~diff([-inf;tmp;inf]))
 112         %ix2=diff(ix)
 113         %MODE(k)= tmp(max(ix2)==ix2)
 114 %        Q0500(k) = flix(tmp,R.N(k)/2   + 0.5);
 115 %        Q0250(k) = flix(tmp,R.N(k)/4   + 0.5);
 116 %        Q0750(k) = flix(tmp,R.N(k)*3/4 + 0.5);
 117 %end;
 118 %R.MEDIAN       = Q0500;
 119 %R.Quartiles    = [Q0250; Q0750];
 120
 121 %R.Skewness.Fisher = (R.CM3)./(R.STD.^3);       %%% same as R.SKEWNESS
 122
 123 %R.Skewness.Pearson_Mode   = (R.MEAN-R.MODE)./R.STD;
 124 %R.Skewness.Pearson_coeff1 = (3*R.MEAN-R.MODE)./R.STD;
 125 %R.Skewness.Pearson_coeff2 = (3*R.MEAN-R.MEDIAN)./R.STD;
 126 %R.Skewness.Bowley = (Q0750+Q0250 - 2*Q0500)./(Q0750-Q0250); % quartile skewness coefficient
 127
 128 R.CM2   = R.SSQ0./n1;
 129 szi = size(i); szm = [size(R.MEAN),1];
 130 i       = i - repmat(R.MEAN,szi./szm(1:length(szi)));
 131 R.CM3   = sumskipnan(i.^3,DIM)./n1;
 132 R.CM4   = sumskipnan(i.^4,DIM)./n1;
 133 %R.CM5  = sumskipnan(i.^5,DIM)./n1;
 134
 135 R.SKEWNESS = R.CM3./(R.STD.^3);
 136 R.KURTOSIS = R.CM4./(R.VAR.^2)-3;
 137 [R.MAD,N]  = sumskipnan(abs(i),DIM);    % mean absolute deviation
 138 R.MAD = R.MAD./n1;
 139
 140 R.datatype = 'STAT Level 3';
 141
 142 tmp = version;
 143 if 0, %str2num(tmp(1))*1000+str2num(tmp(3))*100+str2num(tmp(5:6))<2136,
 144         % ###obsolete: was needed for Octave version < 2.1.36
 145         if strcmp(fun(1:2),'CM')
 146                 oo = str2double(fun(3:length(fun)));
 147                 varargout  = sumskipnan(i.^oo,DIM)./n1;
 148         elseif isempty(fun)
 149                 varargout  = R;
 150         else
 151                 varargout  = getfield(R,upper(fun));
 152         end;
 153 else
 154         if iscell(fun),
 155                 for k=1:length(fun),
 156                         if strcmp(fun{k}(1:2),'CM')
 157                                 oo = str2double(fun{k}(3:length(fun{k})));
 158                                 varargout{k}  = sumskipnan(i.^oo,DIM)./n1;
 159                         else
 160                                 varargout{k}  = getfield(R,upper(fun{k}));
 161                         end;
 162                 end;
 163         elseif ischar(fun),
 164                 if strcmp(fun(1:2),'CM')
 165                         oo = str2double(fun(3:length(fun)));
 166                         varargout{1}  = sumskipnan(i.^oo,DIM)./n1;
 167                 else
 168                         varargout{1}  = getfield(R,upper(fun));
 169                 end;
 170         else
 171                 varargout{1} = R;
 172         end;
 173 end;