www.gusucode.com > 时间序列分析工具箱 - tsa源码程序 > tsa/hist2res.m
function [R]=hist2res(H,fun) % Evaluates Histogram data % [R]=hist2res(H) % % [y]=hist2res(H,fun) % estimates fun-statistic % % fun 'mean' mean % 'std' standard deviation % 'var' variance % 'sem' standard error of the mean % 'rms' root mean square % 'meansq' mean of squares % 'sum' sum % 'sumsq' sum of squares % 'CM#' central moment of order # % 'skewness' skewness % 'kurtosis' excess coefficient (Fisher kurtosis) % % see also: NaN/statistic % % REFERENCES: % [1] C.L. Nikias and A.P. Petropulu "Higher-Order Spectra Analysis" Prentice Hall, 1993. % [2] C.E. Shannon and W. Weaver "The mathematical theory of communication" University of Illinois Press, Urbana 1949 (reprint 1963). % [3] http://www.itl.nist.gov/ % [4] http://mathworld.wolfram.com/ % $Id$ % Copyright (C) 1996-2004 by Alois Schloegl <a.schloegl@ieee.org> % This is part of the BIOSIG-toolbox http://biosig.sf.net/ % 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 2 % 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, write to the Free Software % Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111-1307, USA. if ~strcmp(H.datatype,'HISTOGRAM') fprintf(2,'ERROR: arg1 is not a histogram\n'); return; end; if nargin<2, fun=[]; end; sz = size(H.H)./size(H.X); R.N = sumskipnan(H.H,1); R.SUM = sumskipnan(H.H.*repmat(H.X,sz),1); R.SSQ = sumskipnan(H.H.*repmat(H.X.*H.X,sz),1); %R.S3P = sumskipnan(H.H.*repmat(H.X.^3,sz),1); % sum of 3rd power R.S4P = sumskipnan(H.H.*repmat(H.X.^4,sz),1); % sum of 4th power %R.S5P = sumskipnan(H.H.*repmat(H.X.^5,sz),1); % sum of 5th power 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_estimation, 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; x = repmat(H.X,sz) - repmat(R.MEAN,size(H.X,1),1); R.CM3 = sumskipnan(H.H.*(x.^3),1)./n1; R.CM4 = sumskipnan(H.H.*(x.^4),1)./n1; %R.CM5 = sumskipnan(H.H.*(x.^5),1)./n1; R.SKEWNESS = R.CM3./(R.STD.^3); R.KURTOSIS = R.CM4./(R.VAR.^2)-3; R.MAD = sumskipnan(H.H.*abs(x),1)./R.N; % mean absolute deviation tmp = diff(H.X); R.QUANT = min(tmp); H.PDF = H.H./H.N(ones(size(H.H,1),1),:); R.ENTROPY = -sumskipnan(H.PDF.*log2(H.PDF),1); if ~isempty(fun), fun = upper(fun); if strncmp(fun,'CM',2) oo = str2double(fun(3:length(fun))); R = sumskipnan(H.PDF.*(x.^oo),1); else R = getfield(R,fun); end; end;