function [Theta, Bias, BcTheta, BcThetash] = bsccrie(X,Y,nb);
% bsccrie model
%
% 1. Function name
% ----------------
% bs  : Bootstrapped
% ccr : Charnes, Cooper, and Rhodes
% I   : Input-Oriented
% E   : Envelopment Form
%
% 2. Input arguments
% ------------------
% X : (n x m) input matrix
% Y : (n x s) output matrix
% nb : number of bootstrap replications
%
% 3. Output arguments
% -------------------
% Theta     : (n x 1)  efficiency scores
% Thetash   : (n x nb) bootstrap estimate of Theta
% Bias      : (n x 1)  Bias
% BcTheta   : (n x 1)  Bias-corrected Theta
% BcThetash : (n x nb) Bias-corrected bootstrap estimate of Theta
%  where, n  : number of DMUs   
%         nb : number of bootstrap replications
%
% 4. Example
% ----------
% load Rhodes.dat;
% x = Rhodes(:,1:5);
% y = Rhodes(:,6:8);
% [theta,bias,bctheta,bcthetash] = bsccrie(x,y,1000);
%

% 5. Reference
% ------------
%
% 6. Version
% ----------

tic;
[n,m] = size(X);
[n,s] = size(Y);

Theta = zeros(n,1);
Thetash = zeros(n,nb);
Bias = zeros(n,1);
BcTheta = zeros(n,1);
BcThetash = zeros(n,nb);

XT=X';
YT=Y';

for k = 1:n
	%disp(int2str(k));
	c=[ones(1), zeros(1,n), zeros(1,m) , zeros(1,s)]';
	A=[ XT(:,k)    -XT       -eye(m)      zeros(m,s) ;
	    zeros(s,1)  YT        zeros(s,m) -eye(s)     ];
	b=[ zeros(m,1);
	    YT(:,k)   ];
	csense=repmat('E',1,m+s);
	l=[];
	u=[];
	[x,y,d1,d2,obj,solstat] = LMsolvem(A,b,c,csense,l,u);
	Theta(k) = x(1);
end

mt = mean(Theta);
stdt = std(Theta);
iqrt = iqr(Theta);
h = 0.9*(n.^(-0.2))*min(stdt,iqrt);

% start the boostrap iteration
for bs = 1:nb
	seed = bs;
	rand('state',seed);
	if mod(bs,20) == 0   % control the display of the iteration index
		disp(int2str(bs));
	end
	idx = unidrnd(n,n,1);
	betas = Theta(idx);
	rand('state',seed);
	eps = normrnd(0,1,n,1);
	thetast = zeros(n,1);
	for k = 1:n
		if betas(k) + h * eps(k) <= 1
			thetast(k) = betas(k) + h * eps(k);
		else 
			thetast(k) = 2 - betas(k) - h * eps(k);
		end
	end
	mbetas = mean(betas);
	vtheta = var(Theta);
	thetas = mbetas + (1 ./ (1 + h.^2 ./ vtheta).^0.5) .*(thetast - mbetas);
	XS = repmat((Theta ./ thetas),1,m) .* X;
	XST = XS';
	
	for k = 1:n
		%disp(int2str(k));
		c=[ones(1), zeros(1,n), zeros(1,m) , zeros(1,s)]';
		A=[ XT(:,k)    -XST      -eye(m)      zeros(m,s) ;
		    zeros(s,1)  YT        zeros(s,m) -eye(s)     ];
		b=[ zeros(m,1);
		    YT(:,k)   ];
		csense=repmat('E',1,m+s);
		l=[];
		u=[];
		[x,y,d1,d2,obj,solstat] = LMsolvem(A,b,c,csense,l,u);
		Thetash(k,bs) = x(1);
	end
end

ThetashT = Thetash';
mThetashT = mean(ThetashT);
mThetash = mThetashT';
Bias = mThetash - Theta;
BcTheta = Theta - Bias;

for k = 1:n
	for bs = 1:nb
		BcThetash(k,bs) = Thetash(k,bs) - 2*Bias(k);
	end
end

toc;