1 function [varargout]=statistic(i,DIM,fun)
2 % STATISTIC estimates various statistics at once.
5 % calculates all statistic (see list of fun) in dimension DIM
6 % R is a struct with all statistics
9 % estimate of fun on dimension DIM
10 % y gives the statistic of fun
15 % N: STATS of N-th dimension
16 % default or []: first DIMENSION, with more than 1 element
19 % 'std' standard deviation
21 % 'sem' standard error of the mean
22 % 'rms' root mean square
23 % 'meansq' mean of squares
25 % 'sumsq' sum of squares
26 % 'CM#' central moment of order #
28 % 'kurtosis' excess coefficient (Fisher kurtosis)
29 % 'mad' mean absolute deviation
32 % - can deal with NaN's (missing values)
33 % - dimension argument
34 % - compatible to Matlab and Octave
36 % see also: SUMSKIPNAN
39 % [1] http://www.itl.nist.gov/
40 % [2] http://mathworld.wolfram.com/
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/
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.
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.
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/>.
74 DIM = find(size(i)>1,1);
75 if isempty(DIM), DIM=1; end;
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
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
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
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;
102 q = quantile(i, (1:3)/4, DIM);
104 %sz=size(i);sz(DIM)=1;
105 %Q0500=repmat(nan,sz);
110 % tmp = sort(i(:,k));
111 %ix = find(~~diff([-inf;tmp;inf]))
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);
119 %R.Quartiles = [Q0250; Q0750];
121 %R.Skewness.Fisher = (R.CM3)./(R.STD.^3); %%% same as R.SKEWNESS
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
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;
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
140 R.datatype = 'STAT Level 3';
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;
151 varargout = getfield(R,upper(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;
160 varargout{k} = getfield(R,upper(fun{k}));
164 if strcmp(fun(1:2),'CM')
165 oo = str2double(fun(3:length(fun)));
166 varargout{1} = sumskipnan(i.^oo,DIM)./n1;
168 varargout{1} = getfield(R,upper(fun));