%
%
%                        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 = 4/3;
  n2 = 3;

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

  X2 = [theta1*t' theta2*t'];
  V2 = inv(X2'*X2 + 1e-10*eye(2,2))*var;

  obs2 = y2 + error;

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

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

  [results2,chain2,s2chain] = mcmcrun(model,data,params,options);

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

  t1vals_prob2 = chain2(:,1);
  t2vals_prob2 = chain2(:,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);


%  figure(2); clf
%  mcmcplot(chain2,[],results2,'chainpanel');
  

  figure(2)
  scatter(t1vals_prob2,t2vals_prob2,'b')
  box on
  axis([1 3 1 3])
  axis square
  set(gca,'Fontsize',[23]);
  xlabel('\theta_1')
  ylabel('\theta_2')

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

  figure(4)
  nbins = 20;
  h = histnorm(t1vals_prob2(M-M/4:M),nbins,'plot');
  set(gca,'FontSize',[23])
  xlabel('\theta_1')
  ylabel('Frequency')