%  Assume a statistical model of the form
%
%  y = X(b) + e, e~ndraws(0,sigma^2*I)         (*******)
%
%  where y is nx1, beta mx1, and X:R^m-->R^n is nonlinear. 
%
%  In this file we use several sampling methods to sample from the b_hat
%  distribution. 
%
clear all, close all
%--------------------------------------------------------------------------
%  Enter data and all other problem dependent variables. Then fit the data. 
%--------------------------------------------------------------------------
% Input data and other problem dependent functions and parameters.
tvec = [1 3 5 7 9]';
yvec = [0.076 0.258 0.369 0.492 0.559]';
nlin_fn = @expfun;
b_0 = ones(2,1);
nparams = length(b_0);
npts = length(yvec);

% Fit the data using nlinfit.
[b_opt,r,J,C,mse] = nlinfit(tvec,yvec,nlin_fn,b_0);
model_opt = feval(nlin_fn,b_opt,tvec);
figure(1)
  plot(tvec,yvec,'o',tvec,model_opt)
  title('Data and Model Fit')    

%--------------------------------------------------------------------------
% Compute Random draws from distribution for b_1 and b_2 using Monte Carlo 
% simulation. Then compute nonparametric confidence intervals.
%--------------------------------------------------------------------------
ndraws = 10000; 
sig = sqrt(mse);
b_estMC = zeros(nparams,ndraws);
disp('*** Monte Carlo Sampling ***')
for i=1:ndraws
  h = waitbar(i/ndraws);
  % Monte Carlo sampling method
  eMC = sig*randn(npts,1);       
  y = model_opt+eMC;
  b_estMC(:,i) = nlinfit(tvec,y,nlin_fn,b_opt);
end
close(h)

%--------------------------------------------------------------------------
% Compute Random draws from distribution for b_1 and b_2 using bootstrap
% sampling. Then compute nonparametric confidence intervals.
%--------------------------------------------------------------------------
sig = sqrt(mse);
b_estBS = zeros(nparams,ndraws);
disp('*** Bootstrap Sampling ***')
for i=1:ndraws
  h = waitbar(i/ndraws);
  % Bootstrap sampling method
  eBS = randsample(r,npts,1); 
  y = model_opt+eBS;
  b_estBS(:,i) = nlinfit(tvec,y,nlin_fn,b_opt);
end
close(h)

%--------------------------------------------------------------------------
% Compute Random draws from (*******) above using Markov Chain Monte Carlo 
% simulation and bootstrapping. 
%--------------------------------------------------------------------------
bounds = [zeros(1,2);10000*ones(1,2)];
adapt_int = 200;
C = mse*inv(J'*J);
cost_fn = 'expfun_ls';
cost_params.tvec = tvec;
cost_params.yvec = yvec;
disp('*** MCMC Sampling ***')
b_estMCMC = Metropolis(ndraws,b_opt',C,mse,cost_fn,cost_params,bounds,adapt_int);

%--------------------------------------------------------------------------
% Compute confidence intervals using: normal theory, Monte Carlo
% simulation, bootstrapping, and MCMC.
%--------------------------------------------------------------------------
disp('----------Normal Theory       -----------------------------')
ci = nlparci(b_opt,r,'jacobian',J);
for i=1:nparams
  fprintf('b_%d normal theory: 0.025 quantile = %2.5f, mean = %2.5f 0.975 quantile = %2.5f\n',i,ci(i,1),mean([ci(i,1),ci(i,2)]),ci(i,2))
end
disp('----------Monte Carlo Sampling-----------------------------')
sample_plot(b_estMC,b_opt,2)
disp('----------Bootstrap Sampling  -----------------------------')
sample_plot(b_estBS,b_opt,5)
disp('----------MCMC Sampling       -----------------------------')
sample_plot(b_estMCMC',b_opt,8)