octave_packages/tsa-4.2.4/histo2.m

   1 function R = histo2(Y, W)
   2 % HISTO2 calculates histogram for multiple columns with separate bin values
   3 %    for each data column.
   4 %
   5 % R = HISTO2(Y)
   6 % R = HISTO2(Y, W)
   7 %       Y       data
   8 %       W       weight vector containing weights of each sample,
   9 %               number of rows of Y and W must match.
  10 %               default W=[] indicates that each sample is weighted with 1.
  11 %
  12 % R = HISTO(...)
  13 %       R is    a struct with th fields
  14 %       R.X     the bin-values, bin-values are computed separately for each
  15 %               data column, thus R.X is a matrix, each column contains the
  16 %               the bin values of for each data column, unused elements are indicated with NaN.
  17 %               In order to have common bin values, use HISTO3.
  18 %       R.H  is the frequency of occurence of value X
  19 %       R.N  are the number of valid (not NaN) samples (i.e. sum of weights)
  20 %
  21 % more histogram-based results can be obtained by HIST2RES2
  22 %
  23 % see also: HISTO, HISTO2, HISTO3, HISTO4
  24 %
  25 % REFERENCE(S):
  26 %  C.E. Shannon and W. Weaver "The mathematical theory of communication" University of Illinois Press, Urbana 1949 (reprint 1963).
  27
  28 %       $Id: histo2.m 8383 2011-07-16 20:06:59Z schloegl $
  29 %       Copyright (C) 1996-2002,2008,2011 by Alois Schloegl <alois.schloegl@gmail.com>
  30 %       This is part of the TSA-toolbox
  31 %       http://hci.tugraz.at/~schloegl/matlab/tsa/
  32 %
  33 %    This program is free software: you can redistribute it and/or modify
  34 %    it under the terms of the GNU General Public License as published by
  35 %    the Free Software Foundation, either version 3 of the License, or
  36 %    (at your option) any later version.
  37 %
  38 %    This program is distributed in the hope that it will be useful,
  39 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  40 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  41 %    GNU General Public License for more details.
  42 %
  43 %    You should have received a copy of the GNU General Public License
  44 %    along with this program.  If not, see <http://www.gnu.org/licenses/>.
  45
  46
  47 %%%%% check input arguments %%%%%
  48 [yr,yc] = size(Y);
  49 if nargin < 2,
  50         W = [];
  51 end;
  52 if ~isempty(W) && (yr ~= numel(W)),
  53         error('number of rows of Y does not match number of elements in W');
  54 end;
  55
  56 %%%%% identify all possible X's and generate overall Histogram %%%%%
  57 N  = sum(~isnan(Y), 1);
  58 NN = N;
  59 if isempty(W)
  60         sY  = sort(Y,1);
  61 else
  62         [sY, idx] = sort(Y,1);
  63         W = cumsum(W(idx));     %% W becomes cumulative sum
  64 end;
  65 [ix,iy] = find( diff(sY, [], 1) > 0);
  66 nn0 = 0;
  67
  68 for k = 1:yc,
  69         tmp = [ix(iy==k); N(k)];
  70         nn1 = length(tmp);
  71
  72         if isempty(W)
  73                 H(1:nn1,k) = [tmp(1); diff(tmp)];
  74         else
  75                 %%% Note that W is the cumulative sum
  76                 H(1:nn1,k) = [W(tmp(1),k); diff(W(tmp,k))];
  77                 NN(k) = W(N(k), k);
  78         end;
  79         X(1:nn1, k) = sY(tmp, k);
  80
  81         if k==1;
  82                 nn0 = nn1;
  83         elseif nn1 < nn0,
  84                 H (1+nn1:nn0, k) = NaN;
  85                 X (1+nn1:nn0, k) = NaN;
  86         elseif nn1 > nn0,
  87                 H (1+nn0:nn1, 1:k-1) = NaN;
  88                 X (1+nn0:nn1, 1:k-1) = NaN;
  89                 nn0 = nn1;
  90         end;
  91 end;
  92
  93 R.datatype = 'HISTOGRAM';
  94 R.H = H;
  95 R.X = X;
  96 R.N = NN;
  97