%
% 
%                        id_example_nonlinear.m
%
%
%  This is the first model Example 12.13.
%

  clear all

%
% Input the parameters and set up the model.
%

  theta1 = 2;
  theta2 = 2;
  sigma = 0.1;
  var = sigma^2;
  M = 1e+4;
  n1 = -16;
  n2 = 20;
  prior_x = n1:.1:n2;
  prior_y = 1/(n2-n1)*ones(length(prior_x));;

  t = 0:.01:1;
  error = sigma*randn(size(t));
  n = length(t);
  y1 = (theta1 + theta2)*t;

  X1 = [t' t'];
  V1 = inv(X1'*X1 + 1e-10*eye(2,2))*var;

  obs1 = y1 + error;

%
% Set up DRAM parameters and run DRAM.  This requires the codes in MCMC_STAT.
%

  data.xdata = t;
  data.ydata = obs1;
  params = {
     {'theta1',theta1,n1,n2}
     {'theta2',theta2,n1,n2}}
  model.ssfun = @id1ss;
  model.S20 = sigma;
  model.N0 = 1;
  options.qcov = V1;
  options.nsimu = M;
  options.updatesigma = 1;

  [results1,chain1,s2chain] = mcmcrun(model,data,params,options);

%
% Get results and plot chains and joint sample plots.
%

  t1vals_prob1 = chain1(:,1);
  t2vals_prob1 = chain1(:,2);

%  [~,density_t1,t1mesh,~]=kde(t1vals_prob2(M-1e+3:M),2^4,min(t1vals_prob2)-0.1,max(t1vals_prob2)+0.1);
  
  t1mesh = 1:.05:3.2;
  %[density_t1,t1] = ksdensity(t1vals_prob2(M-1e+3:M),t1mesh,'support',[4/3 3]);
  %[density_t1,t1] = ksdensity(t1vals_prob2(M-1e+3:M),t1mesh);

%
% Compute sample correlation
%

  num = 0;
  den1 = 0;  
  den2 = 0;
  t1_bar = sum(t1vals_prob1)/M;
  t2_bar = sum(t2vals_prob1)/M;
%  for i=1:M
%    num = num + (t1vals_prob1(i,1) - t1_bar)*(t2vals_prob1(i,1) - t2_bar); 
%    den1 = den1 + (t1vals_prob1(i,1) - t1_bar)^2;
%    den2 = den2 + (t2vals_prob1(i,1) - t2_bar)^2;
%  end
  num = sum((t1vals_prob1 - t1_bar).*(t2vals_prob1 - t2_bar));
  den1 = sum((t1vals_prob1 - t1_bar).^2);
  den2 = sum((t2vals_prob1 - t2_bar).^2);

  corr = num/sqrt(den1*den2)


  figure(1); clf
  mcmcplot(chain1,[],results1,'chainpanel');

  figure(2)
  scatter(t1vals_prob1,t2vals_prob1,'b')
  box on
  axis([n1 n2 n1 n2])
  axis square
  hold on
  plot([n1 n2],[0 0],'k-',[0 0],[n1 n2],'k-','linewidth',2)
  hold off
  set(gca,'Fontsize',[23]);
  xlabel('\theta_1')
  ylabel('\theta_2')

  figure(3)
  plot(t1vals_prob1,'-b')
  axis([0 M n1 n2])
  set(gca,'Fontsize',[23])
  xlabel('Chain Iteration')
  ylabel('\theta_1')


  nbins = 20;
  figure(4)
  h = histnorm(t1vals_prob1(M-M/4:1:M),nbins,'plot');
  axis([n1-5 n2+5 0 0.04])
  hold on
  plot(prior_x,prior_y,'-r','linewidth',3)
  hold off 
  set(gca,'FontSize',[23])
  xlabel('\theta_1')
  ylabel('Frequency')
  legend('Posterior','Prior','Location','NorthEast')