function [bootstat, bootsam] = bootstrp(nboot,bootfun,varargin)
%BOOTSTRP Bootstrap statistics.
%   BOOTSTRP(NBOOT,BOOTFUN,...) draws NBOOT bootstrap data samples and
%   analyzes them using the function, BOOTFUN. NBOOT must be a
%   positive integer. The third and later arguments are the data;
%   BOOTSTRP passes bootstrap samples of the data to BOOTFUN.
%
%   [BOOTSTAT,BOOTSAM] = BOOTSTRP(...) Each row of BOOTSTAT contains
%   the results of BOOTFUN on one bootstrap sample. If BOOTFUN
%   returns a matrix, then this output is converted to a long
%   vector for storage in BOOTSTAT. BOOTSAM is a matrix of indices
%   into the rows of the extra arguments.
%
%   Examples:
%   To produce a sample of 100 bootstrapped means of random samples
%   taken from the vector Y:
%      M = BOOTSTRP(100, 'mean', Y)
%
%   To produce a sample of 200 bootstrapped coefficient vectors for
%   a regression of the vector Y on the matrix X:
%      B = BOOTSTRP(200, 'regress', Y, X)

%   Reference:
%      Efron, Bradley, & Tibshirani, Robert, J.
%      "An Introduction to the Bootstrap", 
%      Chapman and Hall, New York. 1993.

%   B.A. Jones 9-27-95, ZP You 8-13-98
%   Copyright 1993-2000 The MathWorks, Inc. 
%   $Revision: 2.10 $  $Date: 2000/05/26 17:28:23 $

% Initialize matrix to identify scalar arguments to bootfun.
la = length(varargin);
scalard = zeros(la,1); db = cell(la,1);

% find out the size information in varargin.
n = 1;
for k = 1:la
   [row,col] = size(varargin{k});
   if max(row,col) == 1
      scalard(k) = 1;
   end
   if row == 1 & col ~= 1
      row = col;
      varargin{k} = varargin{k}(:);
   end
   n = max(n,row);
end

% Create index matrix of bootstrap samples.
bootsam = unidrnd(n,n,nboot);

% Get result of bootfun on actual data and find its size. 
thetafit = feval(bootfun,varargin{:});
[ntheta ptheta] = size(thetafit);

% Initialize a matrix to contain the results of all the bootstrap calculations.
bootstat = zeros(nboot,ntheta*ptheta);

% Do bootfun - nboot times.
for bootiter = 1:nboot
   for k = 1:la
      if scalard(k) == 0
         db{k} = varargin{k}(bootsam(:,bootiter),:);
      else
         db{k} = varargin{k};
      end
   end
   tmp = feval(bootfun,db{:});
   bootstat(bootiter,:) = (tmp(:))';
end
