%%
%                            Example7_15_long.m
%
%
% Construct synthetic data and sampling distribution for the damping
% parameter C as illustrated in Example 7.15.  This code also constructs
% the density from 10,000 simulations, which is compared with the sampling
% distribution.  Because this requires repeated optimization, the code
% takes approximately 40 seconds to run.
%
% Required functions: spring_function.m
%                     ksdensity.m
%

  clear all
  close all

%%
% Specify true parameter and measurement variance values.
%

  m = 1;
  k = 20.5;
  c = 1.5;
  sigma = .1;
  var = sigma^2;
  
%%
% Construct synthetic data, the sensitivity matrix, and covariance matrix V.
%

  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;

%%
% Plot the model response, synthetic data, and residuals.  One notes in the residual plot that the
% 2 sigma interval contains approximately 95% of the residuals.
%

  figure(1)
  plot(t,y,'k','linewidth',2)
  hold on
  plot(t,obs,'x','linewidth',1)
  plot(t,0*y,'k','linewidth',2)
  hold off
  set(gca,'Fontsize',[20]);
  legend('Model','Synthetic Data','Location','NorthEast')
  xlabel('Time (s)')
  ylabel('Displacement')

  figure(2)
  plot(t,error,'x',t,0*y,'k','linewidth',2)
  hold on
  plot(t,2*sigma*ones(size(t)),'--r',t,-2*sigma*ones(size(t)),'--r','linewidth',2)
  hold off
  set(gca,'Fontsize',[20]);
  xlabel('Time (s)')
  ylabel('Residuals')

  
%%
% Construct sampling distribution and determine optimal value of C using the routine fminsearch.m
% for this data.  We additionally compute the confidence interval.
%


  C = [1.4:.001:1.6];
  sample_dist = normpdf(C,1.5,sqrt(V));

  modelfun = @(c) spring_function(c,m,k,obs,t);
  [c_opt,fval] = fminsearch(modelfun,c);

  c_opt
  conf_int = [c_opt - 1.96*sqrt(V) c_opt + 1.96*sqrt(V)]

%%
% Construct the density by repeated sampling.  The final KDE is performed
% using ksdensity. The variable count compiles the number of 95% confidence
% intervals that contain the true parameter value.
%
  count = 0;
  for j = 1:10000
    error = sigma*randn(size(t));
    y = 2*exp((-c/den)*t).*cos((num/den)*t);
    obs = y + error;
    modelfun = @(c) spring_function(c,m,k,obs,t);
    [c_opt,fval] = fminsearch(modelfun,c);
    c_lower = c_opt - 1.96*sqrt(V);
    c_upper = c_opt + 1.96*sqrt(V);
    if (c_lower <= c) & (c <= c_upper)
      count = count + 1;
    end
    Cvals(j) = c_opt;
  end
  count

  [density_C,Cmesh] = ksdensity(Cvals);

%%
% Plot the sampling distribution and constructed density computed by generating data
% and optimizing 10,000 times. 
%
  
  ylabel('Residuals')
 
  figure(3)
  plot(Cmesh,density_C,'k--','linewidth',3)
  axis([1.4 1.6 0 25])
  hold on
  plot(C,sample_dist,'linewidth',3)
  set(gca,'Fontsize',[22]);
  xlabel('Optimal C')
  ylabel('Density')
  legend('Constructed Density','Sampling Density','Location','NorthEast')