%
%  Code to construct prediction intervals for the spring model.
%
  clear all

  global m k

%
% Input the parameters and construct synthetic data
%

  m = 1;
  c = 1.5; 
  k = 20.5;
  sigma = .1;
  var = sigma^2;

  num = sqrt(4*m*k - c^2);
  den = 2*m;

  t = 0:.01:5;
  error = sigma*randn(size(t));
  n = length(t);

  y = 2*exp((-c/den)*t).*cos((num/den)*t);
  dydc = -(t/m).*exp((-c/den)*t).*cos((num/den)*t) + (c*t/(m*num)).*exp((-c/den)*t).*sin((num/den)*t);
  V = inv(dydc*dydc')*var;

  obs = y + error;

%
% Set up DRAM parameters
%

  data.xdata = t;
  data.ydata = obs;
  params = {{'c',c,0}}
  n = 501;
  p =  1;
  
  model.ssfun = @springss;
  model.S20 = sigma;
  model.N0 = 1;
  options.nsimu = 20000;
  options.updatesigma = 1;

%
%  Constuct the parameter density for the damping c.
%

  [results,chain,s2chain] = mcmcrun(model,data,params,options);

%
% Construct and plot prediction intervals.
%

  out = mcmcpred(results,chain,s2chain,t',@spring_fun,5000);
  mcmcpredplot_custom(out)
  hold on
  plot(t,obs,'xr')
  plot(t,0*t,'-k','linewidth',2)
  hold off
  set(gca,'Fontsize',[20])
  xlabel('Time')
  ylabel('Displacement')