`                          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)]
```
```c_opt =

1.4512

conf_int =

1.4154    1.4871

```

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);
```
```count =

9510

```

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')
```