% MonteCarloExpFun.m
%
%  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. 
%
%  The least squares solution b_hat is computed for a large number of 
%  samples from (*******) and
%
%  b_hat = arg_min_b ||X(b)-y||^2.
%
%  In this file we use Monte Carlo simulation 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);

%--------------------------------------------------------------------------
% Compute Random draws from (*******) above using Monte Carlo simulation.
% Then plot samples and model fit of the data.
%--------------------------------------------------------------------------
ndraws = 1000; 
sig = sqrt(mse);
b_est = zeros(nparams,ndraws);
for i=1:ndraws
  h = waitbar(i/ndraws);
  y = model_opt+sig*randn(npts,1);
  b_est(:,i) = nlinfit(tvec,y,nlin_fn,b_opt);
  model_fit(:,i) = feval(nlin_fn,b_est(:,i),tvec);
end
close(h)
figure(1)
  subplot(211)
    plot(tvec,yvec,'o',tvec,model_opt)
    title('Data and Model Fit')
  subplot(212)
    plot(tvec,model_fit,'y'), hold on, plot(tvec,model_opt,'r'), hold off
    title('Model Fit and Sample Fit Envelope')
    
%--------------------------------------------------------------------------
% Plot samples and perform various tests for normalcy.
%--------------------------------------------------------------------------
sample_plot(b_est,b_opt,2)
disp('---------------------------------------')
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(1,1),mean(b_est(i,:)),ci(i,2))
end