octave_packages/nan-2.5.5/quantile.m

   1 function Q=quantile(Y,q,DIM,method)
   2 % QUANTILE calculates the quantiles of histograms and sample arrays.
   3 %
   4 %  Q = quantile(Y,q)
   5 %  Q = quantile(Y,q,DIM)
   6 %     returns the q-th quantile along dimension DIM of sample array Y.
   7 %     size(Q) is equal size(Y) except for dimension DIM which is size(Q,DIM)=length(Q)
   8 %
   9 %  Q = quantile(HIS,q)
  10 %     returns the q-th quantile from the histogram HIS.
  11 %     HIS must be a HISTOGRAM struct as defined in HISTO2 or HISTO3.
  12 %     If q is a vector, the each row of Q returns the q(i)-th quantile
  13 %
  14 % see also: HISTO2, HISTO3, PERCENTILE
  15
  16
  17 %       $Id: quantile.m 9601 2012-02-09 14:14:36Z schloegl $
  18 %       Copyright (C) 1996-2003,2005,2006,2007,2009,2011 by Alois Schloegl <alois.schloegl@gmail.com>
  19 %       This function is part of the NaN-toolbox
  20 %       http://pub.ist.ac.at/~schloegl/matlab/NaN/
  21
  22 %    This program is free software; you can redistribute it and/or modify
  23 %    it under the terms of the GNU General Public License as published by
  24 %    the Free Software Foundation; either version 3 of the License, or
  25 %    (at your option) any later version.
  26 %
  27 %    This program is distributed in the hope that it will be useful,
  28 %    but WITHOUT ANY WARRANTY; without even the implied warranty of
  29 %    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  30 %    GNU General Public License for more details.
  31 %
  32 %    You should have received a copy of the GNU General Public License
  33 %    along with this program; If not, see <http://www.gnu.org/licenses/>.
  34
  35 if nargin<3,
  36         DIM = [];
  37 end;
  38 if isempty(DIM),
  39         DIM = find(size(Y)>1,1);
  40         if isempty(DIM), DIM = 1; end;
  41 end;
  42
  43
  44 if nargin<2,
  45         help quantile
  46
  47 else
  48         [q, rix]  = sort(q(:)');        % sort quantile values
  49         [tmp,rix] = sort(rix);  % generate reverse index
  50
  51         SW = isstruct(Y);
  52         if SW, SW = isfield(Y,'datatype'); end;
  53         if SW, SW = strcmp(Y.datatype,'HISTOGRAM'); end;
  54         if SW,
  55                 [yr, yc] = size(Y.H);
  56                 Q = repmat(nan,length(q),yc);
  57                 if ~isfield(Y,'N');
  58                         Y.N = sum(Y.H,1);
  59                 end;
  60
  61                 for k1 = 1:yc,
  62                         tmp = Y.H(:,k1)>0;
  63                         h = full(Y.H(tmp,k1));
  64                         t = Y.X(tmp,min(size(Y.X,2),k1));
  65
  66                         N = Y.N(k1);
  67                         t2(1:2:2*length(t)) = t;
  68                         t2(2:2:2*length(t)) = t;
  69                         x2 = cumsum(h);
  70                         x(1)=0;
  71                         x(2:2:2*length(t)) = x2;
  72                         x(3:2:2*length(t)) = x2(1:end-1);
  73
  74                         % Q(q < 0 | 1 < q,:) = NaN;  % already done at initialization
  75                         Q(q==0,k1) = t2(1);
  76                         Q(q==1,k1) = t2(end);
  77                         n = 1;
  78                         for k2 = find( (0 < q) & (q < 1) )
  79                                 while (q(k2)*N > x(n)),
  80                                         n=n+1;
  81                                 end;
  82
  83                                 if q(k2)*N==x(n)
  84                                         % mean of upper and lower bound
  85                                         Q(k2,k1) = (t2(n)+t2(n+1))/2;
  86                                 else
  87                                         Q(k2,k1) = t2(n);
  88                                 end;
  89                         end;
  90                         Q = Q(rix,:);   % order resulting quantiles according to original input q
  91                 end;
  92
  93
  94         elseif isnumeric(Y),
  95                 sz = size(Y);
  96                 if DIM>length(sz),
  97                         sz = [sz,ones(1,DIM-length(sz))];
  98                 end;
  99
 100                 f  = zeros(1,length(q));
 101                 f( (q < 0) | (1 < q) ) = NaN;
 102                 D1 = prod(sz(1:DIM-1));
 103                 D3 = prod(sz(DIM+1:length(sz)));
 104                 Q  = repmat(nan,[sz(1:DIM-1),length(q),sz(DIM+1:length(sz))]);
 105                 for k = 0:D1-1,
 106                 for l = 0:D3-1,
 107                         xi = k + l * D1*sz(DIM) + 1 ;
 108                         xo = k + l * D1*length(q) + 1;
 109                         t  = Y(xi:D1:xi+D1*sz(DIM)-1);
 110                         t  = t(~isnan(t));
 111                         N  = length(t);
 112
 113                         if (N==0)
 114                                 f(:) = NaN;
 115                         else
 116                                 t  = sort(t);
 117                                 t2(1:2:2*length(t)) = t;
 118                                 t2(2:2:2*length(t)) = t;
 119                                 x = floor((1:2*length(t))/2);
 120                                 %f(q < 0 | 1 < q) = NaN;  % for efficiency its defined outside loop
 121                                 f(q==0) = t2(1);
 122                                 f(q==1) = t2(end);
 123
 124                                 n = 1;
 125                                 for k2 = find( (0 < q) & (q < 1) )
 126                                         while (q(k2)*N > x(n)),
 127                                                 n = n+1;
 128                                         end;
 129
 130                                         if q(k2)*N==x(n)
 131                                                 % mean of upper and lower bound
 132                                                 f(k2) = (t2(n) + t2(n+1))/2;
 133                                         else
 134                                                 f(k2) = t2(n);
 135                                         end;
 136                                 end;
 137                         end;
 138                         Q(xo:D1:xo + D1*length(q) - 1) = f(rix);
 139                 end;
 140                 end;
 141
 142         else
 143                 fprintf(2,'Error QUANTILES: invalid input argument\n');
 144                 return;
 145         end;
 146
 147 end;
 148
 149 %!assert(quantile(1:10,[.2,.5]),[2.5, 5.5])
 150
 151