% NonlinDemo_MonteCarlo.m
clear all
%  The purpose of this example is to demonstrate that with nonlinear models
%  if the noise is i.i.d. Gaussian, the parameter estimates may not be
%  i.i.d. Gaussian. Here we use a model of the form
%
%  y = X(b) + e, e~ndraws(0,sigma^2*I)         (*******)
%
%  where y is 2x1, beta 2x1, and X:R^2-->R^2. We keep things this dimension
%  for visualization purposes.
%
%  The least squares solution b_hat is computed for a large number of 
%  samples from the above linear model:
%
%  b_hat = arg_min_b ||X(b)-y||^2.
%
%  In this file we use Monte Carlo simulation to get a feel for the b_hat
%  distribution. 
%
%--------------------------------------------------------------------------
%  Define the true parameters b and the matrix X in (*******).
%--------------------------------------------------------------------------
clear all, close all
% Input data and other problem dependent functions and parameters.
disp('1. First we input data and regressor variables, define the nonlinear ') 
disp('   function X(b), and then fit the data, which yields an optimal ')
disp('   parameter values b_opt and model fit. Confidence intervals are')
disp('   computed using nlparci. Hit any key to continue.'), pause

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);

% First run the optimzation routine
[b_opt,r,J,CovDelta_b,mse] = nlinfit(tvec,yvec,nlin_fn,b_0);
model_opt = feval(nlin_fn,b_opt,tvec);
ci = nlparci(b_opt,r,'jacobian',J);

%--------------------------------------------------------------------------
% Compute Random draws from (*******) above using Monte Carlo simulation.
%--------------------------------------------------------------------------
disp('2. Next we compute a bunch of random draws from the data distribution')
disp('   given by y=X(b_opt)+e, e~N(0,sigma^2*I), where sigma^2 is estimated')
disp('   using the mse of the fit. The corresponding parameter estimates are')
disp('   random draws from b_hat=argmin_b||y-Xb||^2. Hit any key to continue.'), pause

ndraws = 1000; 
nparams = length(b_opt);
npts = length(yvec);
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.
%--------------------------------------------------------------------------
disp('3. Now we perform a variety of normality tests: (i) compare the ')
disp('   analytic and sample covariances, (ii) plot the data and parameter') 
disp('   samples, and (iii) output normal quantile plots and histograms. ')
disp('   Hit any key to continue.'), pause 

% Plot of samples.
figure(2)
  plot(b_est(1,:),b_est(2,:),'*')
  hold on, plot(b_opt(1),b_opt(2),'r*'), hold off
  title('Plots of estimates. True value in red.')

% Plot the normal quantiles.
figure(3)
  subplot(121)
    qqplot(b_est(1,:))
    title('b_1')
  subplot(122)
    qqplot(b_est(2,:))
    title('b_2')

%--------------------------------------------------------------------------
% Create normalized histograms of parameter samples and make inference.
%--------------------------------------------------------------------------
disp('4. Finally, we compute normalize histograms to make proper probability')
disp('   density functions, perform the chi^2 test for normality, and use the')
disp('   quantile function to compute confidence intervals, which are compared')
disp('   with the normal theory confidence intervals. Hit any key to continue.'), pause

% b_1
figure(4)
[nb_est1,tb1]=nhist_fn(b_est(1,:));
q=quantile(b_est(1,:),[0.025 0.975]);
[h,p]=chi2gof(b_est(1,:));
subplot(121)
      plot(tb1,nb_est1,'ok'), hold on
      jj = min( find( abs(tb1-q(1))==min(abs(tb1-q(1))) ) );
      tb1_025 = tb1(jj);
      plot(tb1_025,0,'*r')
      mb1 = mean(b_est(1,:)); 
      plot(mb1,0,'*m')
      jj = min( find( abs(tb1-q(2))==min(abs(tb1-q(2))) ) );
      tb1_975 = tb1(jj);
      plot(tb1_975,0,'*g'), hold off
      title(['b_1. Chi^2: ',num2str(h),', p = ',num2str(p)])
      legend('distribution','.025 quantile','mean','.975 quantile')
      fprintf('b_1      quantile: 0.025 quantile = %2.5f, mean = %2.5f 0.975 quantile = %2.5f\n',tb1_025,mb1,tb1_975)
      fprintf('b_1 normal theory: 0.025 quantile = %2.5f, mean = %2.5f 0.975 quantile = %2.5f\n',ci(1,1),mean([ci(1,1),ci(1,2)]),ci(1,2))
      disp('---------------------------------------')
  % b_2
  [nb_est2,tb2]=nhist_fn(b_est(2,:));
  q=quantile(b_est(2,:),[0.025 0.975]);
  [h,p]=chi2gof(b_est(2,:));
  subplot(122)
      plot(tb2,nb_est2,'ok'), hold on
      jj = min( find( abs(tb2-q(1))==min(abs(tb2-q(1))) ) );
      tb2_025 = tb2(jj);
      plot(tb2_025,0,'*r')
      mb2 = mean(b_est(2,:));
      plot(mb2,0,'*m')
      jj = min( find( abs(tb2-q(2))==min(abs(tb2-q(2))) ) );
      tb2_975 = tb2(jj);
      plot(tb2_975,0,'*g'), hold off
      title(['b_2. Chi^2: ',num2str(h),', p = ',num2str(p)])
      fprintf('b_2      quantile: 0.025 quantile = %2.5f, mean = %2.5f 0.975 quantile = %2.5f\n',tb2_025,mb2,tb2_975)
      fprintf('b_2 normal theory: 0.025 quantile = %2.5f, mean = %2.5f 0.975 quantile = %2.5f\n',ci(2,1),mean([ci(2,1),ci(2,2)]),ci(2,2))