X-Git-Url: https://git.creatis.insa-lyon.fr/pubgit/?a=blobdiff_plain;ds=sidebyside;f=octave_packages%2Ftsa-4.2.4%2Fy2res.m;fp=octave_packages%2Ftsa-4.2.4%2Fy2res.m;h=216ec61c1ccbbfaec596a93ad94be84bd55c1206;hb=c880e8788dfc484bf23ce13fa2787f2c6bca4863;hp=0000000000000000000000000000000000000000;hpb=1705066eceaaea976f010f669ce8e972f3734b05;p=CreaPhase.git diff --git a/octave_packages/tsa-4.2.4/y2res.m b/octave_packages/tsa-4.2.4/y2res.m new file mode 100644 index 0000000..216ec61 --- /dev/null +++ b/octave_packages/tsa-4.2.4/y2res.m @@ -0,0 +1,121 @@ +function [R]=y2res(Y) +% Y2RES evaluates basic statistics of a data series +% +% R = y2res(y) +% several statistics are estimated from each column of y +% +% OUTPUT: +% R.N number of samples, NaNs are not counted +% R.SUM sum of samples +% R.MEAN mean +% R.STD standard deviation +% R.VAR variance +% R.Max Maximum +% R.Min Minimum +% ... and many more including: +% MEDIAN, Quartiles, Variance, standard error of the mean (SEM), +% Coefficient of Variation, Quantization (QUANT), TRIMEAN, SKEWNESS, +% KURTOSIS, Root-Mean-Square (RMS), ENTROPY +% + +% $Id: y2res.m 5090 2008-06-05 08:12:04Z schloegl $ +% Copyright (C) 1996-2005,2008 by Alois Schloegl +% This is part of the TSA-toolbox +% http://octave.sourceforge.net/ +% http://www.dpmi.tugraz.at/~schloegl/matlab/tsa/ +% +% This program is free software: you can redistribute it and/or modify +% it under the terms of the GNU General Public License as published by +% the Free Software Foundation, either version 3 of the License, or +% (at your option) any later version. +% +% This program is distributed in the hope that it will be useful, +% but WITHOUT ANY WARRANTY; without even the implied warranty of +% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the +% GNU General Public License for more details. +% +% You should have received a copy of the GNU General Public License +% along with this program. If not, see . + + +[R.SUM, R.N, R.SSQ] = sumskipnan(Y,1); +%R.S3P = sumskipnan(Y.^3,1); +R.S4P = sumskipnan(Y.^4,1); +%R.S5P = sumskipnan(Y.^5,1); + +R.MEAN = R.SUM./R.N; +R.MSQ = R.SSQ./R.N; +R.RMS = sqrt(R.MSQ); +R.SSQ0 = R.SSQ-R.SUM.*R.MEAN; % sum square of mean removed + +if 1,%flag_implicit_unbiased_estim, + n1 = max(R.N-1,0); % in case of n=0 and n=1, the (biased) variance, STD and STE are INF +else + n1 = R.N; +end; + +R.VAR = R.SSQ0./n1; % variance (unbiased) +R.STD = sqrt(R.VAR); % standard deviation +R.SEM = sqrt(R.SSQ0./(R.N.*n1)); % standard error of the mean +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 +R.Coefficient_of_variation = R.STD./R.MEAN; + +R.CM2 = R.SSQ0./n1; + +R.Max = max(Y,[],1); +R.Min = min(Y,[],1); + +%R.NormEntropy = log2(sqrt(2*pi*exp(1)))+log2(R.STD); + +Q0500=repmat(nan,1,size(Y,2)); +Q0250=Q0500; +Q0750=Q0500; +%MODE=Q0500; +for k = 1:size(Y,2), + tmp = sort(Y(:,k)); + Q0250(k) = flix(tmp,R.N(k)/4 + 0.75); + Q0500(k) = flix(tmp,R.N(k)/2 + 0.50); + Q0750(k) = flix(tmp,R.N(k)*3/4 + 0.25); + tmp = diff(tmp); + + pdf = diff([0; find(tmp>0); R.N(k)])/R.N(k); % empirical probability distribution + R.ENTROPY(k) = -sumskipnan(pdf.*log(pdf)); + + tmp = tmp(find(tmp)); + q = min(tmp); + qerror = 0; + if isempty(q), + q = NaN; + else + tmp = tmp/q; + qerror = max(abs(tmp-round(tmp))); + end; + R.QUANT(k) = q; + R.Qerror(k) = qerror; +end; +if any(R.Qerror*1e6>R.QUANT) + warning('(Y2RES) Quantization might not be equidistant') +end; + +R.MEDIAN = Q0500; +R.Quartiles = [Q0250; Q0750]; +% R.IQR = H_spread = [Q0750 - Q0250]; +R.TRIMEAN = [Q0250 + 2*Q0500 + Q0750]/4; + +Y = Y - repmat(R.MEAN,size(Y)./size(R.MEAN)); +R.CM3 = sumskipnan(Y.^3,1)./n1; +R.CM4 = sumskipnan(Y.^4,1)./n1; +%R.CM5 = sumskipnan(Y.^5,1)./n1; + +R.SKEWNESS = R.CM3./(R.STD.^3); +R.KURTOSIS = R.CM4./(R.VAR.^2)-3; + +%R.Skewness.Fisher = (R.CM3)./(R.STD.^3); %%% same as R.SKEWNESS + +%R.Skewness.Pearson_Mode = (R.MEAN-R.MODE)./R.STD; +%R.Skewness.Pearson_coeff1 = (3*R.MEAN-R.MODE)./R.STD; +R.Skewness.Pearson_coeff2 = (3*R.MEAN-R.MEDIAN)./R.STD; +R.Skewness.Bowley = (Q0750+Q0250 - 2*Q0500)./(Q0750-Q0250); % quartile skewness coefficient + +R.datatype = 'STAT Level 4'; +