请高手用Matlab程序写一段实现参数Bootstrap方法程序，谢谢！

function f=Gamma_bootstrap(k,gamma,lambda,tao,N,n,B,arfa)
X= gamrnd(k,lambda,1,N);
Y= gamrnd(gamma,tao,1,N);
X_bar=mean(X);
Y_bar=mean(Y);
lambda_bar=k/X_bar;
tao_bar=gamma/Y_bar;
rou_hat=tao_bar/lambda_bar;
R_hat=1-fcdf(gamma/k/rou_hat,2*k,2*gamma);
for i=1:B
XX= gamrnd(k,lambda_bar,1,n);
YY= gamrnd(gamma,tao_bar,1,n);
X_star=mean(XX);
Y_star=mean(YY);
lambda_star(i)=k/X_star;
tao_star(i)=gamma/Y_star;
rou_star(i)=tao_star(i)./lambda_star(i);
R_star(i)=1-fcdf(gamma/k./rou_star(i),2*k,2*gamma);
end
% Lambda=var(lambda_star)
% MSE_Lambda=var(lambda_star)*(n-1)/n+(mean(lambda_star)-lambda)^2
% Tao=var(tao_star)
% MSE_Tao=var(tao_star)*(n-1)/n+(mean(tao_star)-tao)^2
rou=tao/lambda;
real_R=1-fcdf(gamma/k./rou,2*k,2*gamma)
R_star_mean=mean(R_star)
R_star_variance=var(R_star)
R_star_sort=sort(R_star);
R_quantile_left=R_star_sort(round((B+1)*(1-arfa/2)));
R_quantile_right=R_star_sort(round((B+1)*arfa/2));
left_R=2*R_hat-R_quantile_left
right_R=2*R_hat-R_quantile_right
length=right_R-left_R