%% A case study of bootstrap for parameter estimation with Matlab
% Michaelis Mentan Model: V(S)=Vm*S/(Km+S);
clear
% close all
% data =xlsread('MM_data.xlsx');
load MM
data=MM;
% Co=100;
% k=0.03;
% beta0(1)=Co;
% beta0(2)=k;

%% initial guesses
Co=.3;
k=460;
beta0(1)=Co;
beta0(2)=k;
t=data(:,1);
n=length(t);
p=length(beta0);
Cobs=data(:,2);  

%% nlinfit returns beta, residuals, Jacobian (sensitivity coefficient matrix),
%covariance matrix, and mean square error
MM2=@(beta,S) beta(1)*S./(beta(2)+S);
[beta,resids,J,COVB,mse] = nlinfit(t,Cobs,MM2,beta0);
rmse=sqrt(mse);

%% R is the correlation matrix for the betaeters, sigma is the standard deviation vector
[R,sigma]=corrcov(COVB);

%% asymptotic confidence intervals for beta
% 参数置信区间
ci95=nlparci(beta,resids,J); %95% CI,alpha is 0.05 by default
ci90=nlparci(beta,resids,J, 0.1); %90% CI

%% nonlinear regression confidence intervals-- 'on' means simultaneous
%bounds; 'off' is for nonsimultaneous bounds; must use 'curve' for
%regression line, 'observation' for prediction interval
% 回归与预测置信带宽
[Cpred, delta] = nlpredci(MM2,t,beta,resids,J,0.05,'on','curve'); %confidence band for regression line
[Cpred, deltaob] =nlpredci(MM2,t,beta,resids,J,0.05,'on','observation');%prediction band for individual points

%% simultaneous confidence bands for regression line
CBu=Cpred+delta;
CBl=Cpred-delta;

%% simultaneous prediction bands for regression line
PBu=Cpred+deltaob;
PBl=Cpred-deltaob;

%% bootstrap CI for beta
% 两种bootstrap方法：bootstrapping the data和bootstrapping the residuals
nboot=1000; % bootstrap次数
mm=1;%use x1, y1; x2, y2;... bootstrapping the data
% mm=2;%use x1, Ypred1+e1; x2, Ypred2+e2;... bootstrapping the residuals
for j=1:nboot
   r=round(1 + (n-1).*rand(n,1));% index of random integers from 1 to n
   % rand函数creats the n random numbers from a uniform distribution
   % 构造随机随机数以对原始数据有放回取样重组
   for i=1:n % n=21
      if mm==1 % bootstrapping the raw data
          tt(i)=t(r(i));% tt(i) is the time for each bootstrapped datum
          yboot(i)=Cobs(r(i));% yboot(i) is the y value for each bootstrapped datum
      end
      if mm==2 % bootstrappingg the residuals
          tt=t;
          yboot(i)=Cpred(i)+resids(r(i)); %利用bootstrap残差由预测值重构观测值
          if i==n
              yboot=yboot'; %make yboot a column
          end
      end
   end
   % 每bootstrap loop一次，nlinfit参数估计一次
   % collect all estimates of paramete and prediction residuals
   [betab(j,:),rr(j,:),J2,COVB2,mse2]= nlinfit(tt,yboot,MM2,beta0);%betab are the parameters from the bootstrapping
   % betab:1000*2; each column collects corresponding estimates 1000 times
   % rr:1000*21
   ypredb(j,:)=MM2(betab(j,:),t); % records the predicted values based on the bootstrap-based parameters each time
   % ypredb:1000*21
   clear yboot
end

r2=rr(1,:)';%save the residuals from first bootstrapping 提取第一次boostra所得的预测残差
for j=2:nboot
   r2=[r2; rr(j,:)']; %build a matrix "r2" where each column holds residuals for one bootstrap
end
% 残差向量汇与一列，以便绘制残差histogram

% compute the 95% confidence interval(CI) for bootstrapped parameter estimtes
% sort each parameter from lowest to highest
bsort=sort(betab,1); %sorts along columns
% sort(x,1)对x的每一列进行排序，sort(x,2)对x的每一行进行排序，默认ascending order
L=round(0.025*nboot); %L is lower 2.5% cutoff
if L==0;  L=1;  end %if L = 0, switch to L = 1 to avoid MATLAB error
U=round(0.975*nboot); %U is upper 2.5% cutoff
cib(1,1)=bsort(L,1); cib(1,2)=bsort(U,1);%bootstrap 95% CI for first beta
cib(2,1)=bsort(L,2); cib(2,2)=bsort(U,2);%bootstrap 95% CI for second beta
% cib:2*2, similar to ci95
% For y 95% bootstrap confidence bands 回归线置信带
ysort=sort(ypredb,1);
for i=1:n
   ybci(i,1)=ysort(L,i); ybci(i,2)=ysort(U,i);%ybci is a n-by-2 matrix with bootstrap CI for y at each time
end

%% compute bootstrap prediction bands
% Prediction width=sqrt(confidence with^2+(Tcrit*RMSE)^2)
D=rmse*tinv(.975,n-p);
% 先计算回归线置信带宽
% Cpred 不一定落在ybci的中间位置
CIwb(:,1)=Cpred-ybci(:,1); CIwb(:,2)=ybci(:,2)-Cpred; %upper (column 1) and lower (column 2) bootstrap CIwidths
% 预测线置信带宽
PIwb(:,1)=sqrt(CIwb(:,1).^2+D^2); PIwb(:,2)=sqrt(CIwb(:,2).^2+D^2);%upper and lower widths of PI
PIb(:,1)=Cpred-PIwb(:,1); PIb(:,2)=Cpred+PIwb(:,2); %PI values

%% residual histogram for bootstrap residuals
[n1, xout] = hist(r2,6); % r2: 21*1000
figure
hold on
set(gca, 'fontsize',14,'fontweight','bold');
bar(xout, n1) % plots the histogram
xlabel('Observed C - Predicted C','fontsize',16,'fontweight','bold')
ylabel('Frequency','fontsize',16,'fontweight','bold')
%Monte-Carlo

%% plot confidence and prediction bands for regression line
%%plot Cobs, Cpred line, confidence band for regression line
figure
hold on
set(gca, 'fontsize',14,'fontweight','bold');
L4 = ['S (mmol)'];
xlabel(L4,'fontsize',16,'fontweight','bold');
% ylabel('log{itS}_a','fontsize',16,'fontweight','bold');
ylabel('V (mmol/S)','fontsize',16,'fontweight','bold');
h1(1)=plot(t,Cobs,'square', 'Markerfacecolor', 'b');
h1(2) = plot(t,Cpred,'-','LineWidth',2); % predicted values
h1(3) = plot(t,CBu,'--r','LineWidth',2); % confidence upper values
plot(t,CBl,'--r','LineWidth',2); % confidence lower values
h1(4) = plot(t,PBu,'-.','LineWidth',2);
plot(t,PBl,'-.','LineWidth',2);

% plot bootstrap bands
h1(5) = plot(t,ybci(:,1),'--k','LineWidth',2); % bootstrapped y confidence lower values
plot(t,ybci(:,2),'--k','LineWidth',2); % bootstrapped y confidence upper values
h1(6) = plot(t,PIb(:,1),'-g','LineWidth',2);
plot(t,PIb(:,2),'-g','LineWidth',2)
legend(h1,'Cobs','Cpred','asyCB','asyPB','bootCB','bootPB','location','best')