octave_packages/tsa-4.2.4/y2res.m

   1 function [R]=y2res(Y)
   2 % Y2RES evaluates basic statistics of a data series
   3 %
   4 % R = y2res(y)
   5 %       several statistics are estimated from each column of y
   6 %
   7 % OUTPUT:
   8 %   R.N     number of samples, NaNs are not counted
   9 %   R.SUM   sum of samples
  10 %   R.MEAN  mean
  11 %   R.STD   standard deviation
  12 %   R.VAR   variance
  13 %   R.Max   Maximum
  14 %   R.Min   Minimum
  15 %   ...   and many more including:
  16 %       MEDIAN, Quartiles, Variance, standard error of the mean (SEM),
  17 %       Coefficient of Variation, Quantization (QUANT), TRIMEAN, SKEWNESS,
  18 %       KURTOSIS, Root-Mean-Square (RMS), ENTROPY
  19 %
  20
  21 %       $Id: y2res.m 5090 2008-06-05 08:12:04Z schloegl $
  22 %       Copyright (C) 1996-2005,2008 by Alois Schloegl <a.schloegl@ieee.org>
  23 %       This is part of the TSA-toolbox
  24 %       http://octave.sourceforge.net/
  25 %       http://www.dpmi.tugraz.at/~schloegl/matlab/tsa/
  26 %
  27 %    This program is free software: you can redistribute it and/or modify
  28 %    it under the terms of the GNU General Public License as published by
  29 %    the Free Software Foundation, either version 3 of the License, or
  30 %    (at your option) any later version.
  31 %
  32 %    This program is distributed in the hope that it will be useful,
  33 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  34 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  35 %    GNU General Public License for more details.
  36 %
  37 %    You should have received a copy of the GNU General Public License
  38 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  39
  40
  41 [R.SUM, R.N, R.SSQ] = sumskipnan(Y,1);
  42 %R.S3P = sumskipnan(Y.^3,1);
  43 R.S4P = sumskipnan(Y.^4,1);
  44 %R.S5P = sumskipnan(Y.^5,1);
  45
  46 R.MEAN  = R.SUM./R.N;
  47 R.MSQ   = R.SSQ./R.N;
  48 R.RMS   = sqrt(R.MSQ);
  49 R.SSQ0  = R.SSQ-R.SUM.*R.MEAN;          % sum square of mean removed
  50
  51 if 1,%flag_implicit_unbiased_estim,
  52     n1  = max(R.N-1,0);                 % in case of n=0 and n=1, the (biased) variance, STD and STE are INF
  53 else
  54     n1  = R.N;
  55 end;
  56
  57 R.VAR   = R.SSQ0./n1;                   % variance (unbiased)
  58 R.STD   = sqrt(R.VAR);                  % standard deviation
  59 R.SEM   = sqrt(R.SSQ0./(R.N.*n1));      % standard error of the mean
  60 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
  61 R.Coefficient_of_variation = R.STD./R.MEAN;
  62
  63 R.CM2   = R.SSQ0./n1;
  64
  65 R.Max   = max(Y,[],1);
  66 R.Min   = min(Y,[],1);
  67
  68 %R.NormEntropy = log2(sqrt(2*pi*exp(1)))+log2(R.STD);
  69
  70 Q0500=repmat(nan,1,size(Y,2));
  71 Q0250=Q0500;
  72 Q0750=Q0500;
  73 %MODE=Q0500;
  74 for k = 1:size(Y,2),
  75         tmp = sort(Y(:,k));
  76         Q0250(k) = flix(tmp,R.N(k)/4   + 0.75);
  77         Q0500(k) = flix(tmp,R.N(k)/2   + 0.50);
  78         Q0750(k) = flix(tmp,R.N(k)*3/4 + 0.25);
  79         tmp = diff(tmp);
  80
  81         pdf   = diff([0; find(tmp>0); R.N(k)])/R.N(k); % empirical probability distribution
  82         R.ENTROPY(k) = -sumskipnan(pdf.*log(pdf));
  83
  84         tmp = tmp(find(tmp));
  85         q   = min(tmp);
  86         qerror = 0;
  87         if isempty(q),
  88                 q = NaN;
  89         else
  90                 tmp = tmp/q;
  91                 qerror = max(abs(tmp-round(tmp)));
  92         end;
  93         R.QUANT(k) = q;
  94         R.Qerror(k) = qerror;
  95 end;
  96 if any(R.Qerror*1e6>R.QUANT)
  97         warning('(Y2RES) Quantization might not be equidistant')
  98 end;
  99
 100 R.MEDIAN        = Q0500;
 101 R.Quartiles     = [Q0250; Q0750];
 102 % R.IQR = H_spread      = [Q0750 - Q0250];
 103 R.TRIMEAN       = [Q0250 + 2*Q0500 + Q0750]/4;
 104
 105 Y       = Y - repmat(R.MEAN,size(Y)./size(R.MEAN));
 106 R.CM3   = sumskipnan(Y.^3,1)./n1;
 107 R.CM4   = sumskipnan(Y.^4,1)./n1;
 108 %R.CM5  = sumskipnan(Y.^5,1)./n1;
 109
 110 R.SKEWNESS = R.CM3./(R.STD.^3);
 111 R.KURTOSIS = R.CM4./(R.VAR.^2)-3;
 112
 113 %R.Skewness.Fisher = (R.CM3)./(R.STD.^3);       %%% same as R.SKEWNESS
 114
 115 %R.Skewness.Pearson_Mode   = (R.MEAN-R.MODE)./R.STD;
 116 %R.Skewness.Pearson_coeff1 = (3*R.MEAN-R.MODE)./R.STD;
 117 R.Skewness.Pearson_coeff2 = (3*R.MEAN-R.MEDIAN)./R.STD;
 118 R.Skewness.Bowley = (Q0750+Q0250 - 2*Q0500)./(Q0750-Q0250); % quartile skewness coefficient
 119
 120 R.datatype = 'STAT Level 4';
 121