% script for Exercise 3.4 (text 2.1 (a,d,g) and (j): Find 95% conf.int. for parms
clear
% get the data
x = [6 10.1333 14.2667 18.4 22.5333 26.6667]';
t = [3.4935 4.2853 5.1374 5.8181 6.8632 8.1841]';
m = length(x);
n = 2;
% the model: t_i = t_0 + s_2*x_i, where t_i,x_i incoming time distances
% and intercept time t_0 and slowness s_2 of 2nd layer are unknown.
% (a): Find least squares solution for the model parameters t_0, s_2.
G = [ones(m,1),x];
soln = G\t;
resid = t - G*soln;
disp('Intercept time: ')
disp(soln(1));
disp('Slowness in second layer: ')
disp(soln(2));
% the graph
hold on
xnodes = 5:.1:27;
tnodes = soln(1)+xnodes*soln(2);
plot(xnodes,tnodes);
plot(x,t,'o');
title('Least Squares Solution and Data Plot');
figure
plot(x,resid,'o');
title('Residuals Plot for Least Squares Solution');
% (d) Assuming data error is normally distributed with std dev sigma = 0.1
% calculate p-value for this model.
sigma = 0.1;
chiobs = (norm(resid)/sigma).^2
disp('p-value for this model: ');
disp(1-chis_cdf(chiobs,m-n));
% (g): Use IRLS to find 1-norm estimates for t_0 and s_2, plot as above and compare
soln1 = irls(G,t,1e-6,1e-6,1,200);
resid1 = t - G*soln1;
disp('Intercept time: ')
disp(soln(1));
disp('Slowness in second layer: ')
disp(soln(2));
% the graph
figure
hold on
xnodes = 5:.1:27;
tnodes = soln1(1)+xnodes*soln1(2);
plot(xnodes,tnodes);
plot(x,t,'o');
title('1-Norm Solution and Data Plot');
figure
plot(x,resid1,'o');
title('Residuals Plot for 1-Norm Solution');
% (j): Find 95% condfidence interval for parameters t_0 and s_2.
% here we use the least squares results from (a) -- the std dev of
% the error in data is known, so we can compute the std dev of the
% parameters as the sqrt of diagonal of the covariance matrix
cov = sigma^2*inv(G'*G)
% so we obtain confidence intervals with lower/upper bounds
lower = soln-1.96*sqrt(diag(cov))
upper = soln+1.96*sqrt(diag(cov))


