octave_packages/nan-2.5.5/hist2res.m

   1 function [R]=hist2res(H,fun)
   2 % Evaluates Histogram data
   3 % [R]=hist2res(H)
   4 %
   5 % [y]=hist2res(H,fun)
   6 %       estimates fun-statistic
   7 %
   8 % fun   'mean'  mean
   9 %       'std'   standard deviation
  10 %       'var'   variance
  11 %       'sem'   standard error of the mean
  12 %       'rms'   root mean square
  13 %       'meansq' mean of squares
  14 %       'sum'   sum
  15 %       'sumsq' sum of squares
  16 %       'CM#'   central moment of order #
  17 %       'skewness' skewness
  18 %       'kurtosis' excess coefficient (Fisher kurtosis)
  19 %
  20 % see also: NaN/statistic
  21 %
  22 % REFERENCES:
  23 % [1] C.L. Nikias and A.P. Petropulu "Higher-Order Spectra Analysis" Prentice Hall, 1993.
  24 % [2] C.E. Shannon and W. Weaver "The mathematical theory of communication" University of Illinois Press, Urbana 1949 (reprint 1963).
  25 % [3] http://www.itl.nist.gov/
  26 % [4] http://mathworld.wolfram.com/
  27
  28 % This program is free software; you can redistribute it and/or
  29 % modify it under the terms of the GNU General Public License
  30 % as published by the Free Software Foundation; either version 2
  31 % of the License, or (at your option) any later version.
  32 %
  33 % This program is distributed in the hope that it will be useful,
  34 % but WITHOUT ANY WARRANTY; without even the implied warranty of
  35 % MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  36 % GNU General Public License for more details.
  37 %
  38 % You should have received a copy of the GNU General Public License
  39 % along with this program; if not, write to the Free Software
  40 % Foundation, Inc., 51 Franklin Street - Fifth Floor, Boston, MA 02110-1301, USA.
  41
  42 %       $Id: hist2res.m 9387 2011-12-15 10:42:14Z schloegl $
  43 %       Copyright (c) 1996-2002,2006 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
  48 if strcmp(H.datatype,'HISTOGRAM'),
  49
  50 elseif strcmp(H.datatype,'qc:histo')
  51         HDR = H;
  52         if isfield(H,'THRESHOLD'),
  53                 TH  = H.THRESHOLD;
  54         else
  55                 TH = repmat([-inf,inf],HDR.NS,1);
  56         end;
  57         HIS = H.HIS;
  58
  59         % remove overflowing samples
  60         HIS.N = sumskipnan(HIS.H);
  61         for k = 1:size(HIS.H,2);
  62                 t = HIS.X(:,min(k,size(HIS.X,2)));
  63                 HIS.H(xor(t<=min(TH(k,:)), t>=max(TH(k,:))),k) = 0;
  64         end;
  65         Nnew = sumskipnan(HIS.H);
  66         R.ratio_lost = 1-Nnew./HIS.N;
  67         HIS.N = Nnew;
  68
  69         % scale into physical values
  70         if H.FLAG.UCAL,
  71                 %t = HIS.X;
  72                 %for k=1:length(HDR.InChanSelect),
  73                 %       HIS.X(:,k) = t(:,min(size(t,2),k))*HDR.Calib(k+1,k)+HDR.Calib(1,k);
  74                 %end;
  75                 HIS.X = [ones(size(HIS.X,1),1),repmat(HIS.X,1,size(HIS.H,2)./size(HIS.X,2))]*H.Calib;
  76         end;
  77         H = HIS;
  78 else
  79         fprintf(2,'ERROR: arg1 is not a histogram\n');
  80         return;
  81 end;
  82 if nargin<2, fun=[]; end;
  83
  84 global FLAG_implicit_unbiased_estimation;
  85 %%% check whether FLAG was already defined
  86 if ~exist('FLAG_implicit_unbiased_estimation','var'),
  87         FLAG_implicit_unbiased_estimation=[];
  88 end;
  89 %%% set DEFAULT value of FLAG
  90 if isempty(FLAG_implicit_unbiased_estimation),
  91         FLAG_implicit_unbiased_estimation=logical(1);
  92 end;
  93
  94 sz      = size(H.H)./size(H.X);
  95 R.N     = sumskipnan(H.H,1);
  96 R.SUM   = sumskipnan(H.H.*repmat(H.X,sz),1);
  97 R.SSQ   = sumskipnan(H.H.*repmat(H.X.*H.X,sz),1);
  98 %R.S3P  = sumskipnan(H.H.*repmat(H.X.^3,sz),1); % sum of 3rd power
  99 R.S4P   = sumskipnan(H.H.*repmat(H.X.^4,sz),1); % sum of 4th power
 100 %R.S5P  = sumskipnan(H.H.*repmat(H.X.^5,sz),1); % sum of 5th power
 101
 102 R.MEAN  = R.SUM./R.N;
 103 R.MSQ   = R.SSQ./R.N;
 104 R.RMS   = sqrt(R.MSQ);
 105 R.SSQ0  = R.SSQ-R.SUM.*R.MEAN;          % sum square of mean removed
 106
 107 if FLAG_implicit_unbiased_estimation,
 108     n1  = max(R.N-1,0);                 % in case of n=0 and n=1, the (biased) variance, STD and STE are INF
 109 else
 110     n1  = R.N;
 111 end;
 112
 113 R.VAR   = R.SSQ0./n1;                   % variance (unbiased)
 114 R.STD   = sqrt(R.VAR);                  % standard deviation
 115 R.SEM   = sqrt(R.SSQ0./(R.N.*n1));      % standard error of the mean
 116 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
 117 R.Coefficient_of_variation = R.STD./R.MEAN;
 118
 119 R.CM2   = R.SSQ0./n1;
 120 x       = repmat(H.X,sz) - repmat(R.MEAN,size(H.X,1),1);
 121 R.CM3   = sumskipnan(H.H.*(x.^3),1)./n1;
 122 R.CM4   = sumskipnan(H.H.*(x.^4),1)./n1;
 123 %R.CM5  = sumskipnan(H.H.*(x.^5),1)./n1;
 124
 125 R.SKEWNESS = R.CM3./(R.STD.^3);
 126 R.KURTOSIS = R.CM4./(R.VAR.^2)-3;
 127 R.MAD = sumskipnan(H.H.*abs(x),1)./R.N; % mean absolute deviation
 128
 129 H.PDF = H.H./H.N(ones(size(H.H,1),1),:);
 130 status=warning('off');
 131 R.ENTROPY = -sumskipnan(H.PDF.*log2(H.PDF),1);
 132 warning(status);
 133 R.QUANT = repmat(min(diff(H.X,[],1)),1,size(H.H,2)/size(H.X,2));
 134 R.MAX = max(H.X);
 135 R.MIN = min(H.X);
 136 R.RANGE = R.MAX-R.MIN;
 137
 138 if ~isempty(fun),
 139         fun=upper(fun);
 140         if strncmp(fun,'CM',2)
 141                 oo = str2double(fun(3:length(fun)));
 142                 R = sumskipnan(H.PDF.*(x.^oo),1);
 143         else
 144                 R = getfield(R,fun);
 145         end;
 146 end;
 147