bstrag(b,fun,c,x,varargin)

function[res] = bstrag(b,fun,c,x,varargin)
% BSTRAG     Bootstrap a statistic in a grouped sample 
% INPUTS   : b      - number of bootstrap samples 
%            fun    - name of function computing desired statistic(s), with
%                     output  placed in p*q array  or structure. Input argu-
%                     ments are  passed to  fun  in the order in which they
%                     appear in bstrag(..);  thus, fun must have c and x as
%                     arguments 1 and 2.
%            c      - n*1 grouping vector,  with c(i)  identifying the group
%                     of ith row of x(,y,z,..). Any non-NaN values are accep-
%                     table; groups are identified by  distinct values in c.
%                     RESAMPLING IS CONDUCTED SEPARATELY FOR EACH GROUP;  if
%                     group m includes  k rows of x(,y,z,..), any  bootstrap
%                     sample will include k randomly selected rows belonging
%                     to it.
%            x      - n*k data matrix, with observations in rows
%            y,z,.. - (optional) additional  vectors/matrices, referenced by
%                     fun, including  scalars. Input arguments are passed to
%                     fun in the  order  they appear in bstrag(..). Non-sca-
%                     lar inputs must have the same number of rows as x.
% OUTPUTS  : res    - statistic values in original and constructed samples
%                     p*q*(b+1) array if fun outputs an array 
%                     1*(b+1) structure array if fun outputs a structure 
%                     First element of res contains value of fun, computed 
%                     in the actual sample.
% EXAMPLE  : See BSTRAG_DEMO 
% SEE ALSO : BSTRAP, BSTATS; JACKKNIFE, BOOTSTRP (Statistics Toolbox)
% AUTHOR   : Dimitri Shvorob, [e-mail], 3/25/07

if nargin < 1
   error('Input argument "b" is undefined')
end
if nargin < 2
   error('Input argument "fun" is undefined') 
end    
if nargin < 3    
   error('Input argument "c" is undefined')
end
if nargin < 4    
   error('Input argument "x" is undefined')
end
if ~isnumeric(b)
   error('Input argument "b" must be numeric')
end  
if b ~= floor(b) 
   error('Input argument "b" must be an integer')
end
if b <= 0
   error('Input argument "b" must be positive')
end
if ~ischar(fun)
   error('Input argument "fun" must be a string')
end  
if ~isnumeric(x)
   error('Input argument "x" must be numeric')
end  
if ~isnumeric(c)
   error('Input argument "c" must be numeric')
end  
k = nargin - 2;
Z = cell(k,1);
N = nan(k,1);
evalString = 'resi = feval(fun';
for j = 1:k
    switch j
        case 1,    Z{j} = c;
        case 2,    Z{j} = x;
        otherwise, Z{j} = varargin{j-2};
    end        
    evalString = [evalString ',Z' num2str(j) 'boot'];
    if ~isnumeric(Z{j})
       error('Input arguments other than "fun" must be numeric')
    end   
    N(j) = length(Z{j});
end 
evalString = [evalString ');'];
n = unique(N(N > 1));
if ~isscalar(n)
   error('Non-scalar input arguments must have the same number of rows')
end   
u = unique(c);
m = length(u);
M = nan(m,1);
for j = 1:m
    ij   = c == u(j);
    M(j) = sum(ij);
    eval(['ic' num2str(j) ' = find(ij);'])
end
for i = 0:b
    if i == 0
       iboot = (1:n)'; %#ok
    else
       iboot = [];
       for j = 1:m
           nj = M(j);
           ibootj = ceil(nj*rand(nj,1)); %#ok
           ibootj = eval(['ic' num2str(j) '(ibootj);']);
           iboot  = [iboot; ibootj];
       end 
    end
    for j = 1:k
        if N(j) == 1
           eval(['Z' num2str(j) 'boot = Z{j};'])
        else
           eval(['Z' num2str(j) 'boot = Z{j}(iboot,:);']) 
        end 
    end 
    try
       eval(evalString)
    catch
       error('Function "fun" could not be evaluated')  
    end
    if i == 0
       if isstruct(resi) 
          outputStructure = true;
          f = fieldnames(resi);
          r = length(f);
       else
          outputStructure = false; 
          [n1,n2] = size(resi);
          res = nan(n1,n2,n+1);
       end    
    end
    if outputStructure
       for j = 1:r
           eval(['res(i+1).' f{j} ' = resi.' f{j} ';']);
       end
    else
       res(:,:,i+1) = resi;
    end
end