% % Example11_7.m and Example11_14.m % % % Construct a quadratic model for the height. Compute the parameter estimates, % covariance matrix 2-sigma intervals and 95% confidences intervals % clear all % close all load HeightAll.dat for j=1:26 ind = 9; boy(j,:,:) = HeightAll(1+ind*(j-1):ind*j,:); end % % Consider boys 6(linear), 9(cubic), 11, 14***, 15***, 24(linear) % m = 15; n = 9; p = 3; figure(1) plot(boy(m,:,1),boy(m,:,2),'x') xt = boy(m,:,1); yt = boy(m,:,2); x = xt'; y = yt'; % x = [-1 -.7479 -.463 -.1643 -.0027 .2466 .5562 .7781 .9945]'; x = [1 2 3 4 5 6 7 8 9]'; % y = [145.8 147.3 148.7 149.78 150.22 152.5 154.8 156.4 158.7]'; % Boy 5 y = [145.5 146.2 148.2 150.3 152 152.3 154.3 156.2 156.8]'; % Boy 18 Pretty iid ... % y = [144.5 146 147.4 149.2 150.8 152.5 155 156.8 158.8]'; % Boy 23 Very good fit but not very iid X = [ones(n,1) x x.^2]; % % Compute the parameter estimates theta, covariance matrix V, error variance % estimate S, and matrix H % theta = inv(X'*X)*(X')*y Model = X*theta; R = y - X*theta; S2 = (1/(n-p))*(R'*R); S = sqrt(S2); V = S2*inv(X'*X); F = inv(X'*X); Hmat = X*inv(X'*X)*(X'); % % Compute 2-sigma intervals % int02s = [theta(1) - 2*sqrt(V(1,1)) theta(1) + 2*sqrt(V(1,1))] int12s = [theta(2) - 2*sqrt(V(2,2)) theta(2) + 2*sqrt(V(2,2))] int22s = [theta(3) - 2*sqrt(V(3,3)) theta(3) + 2*sqrt(V(3,3))] % % Compute 95% confidence intervals using a t-value constructed using tinv.m % alpha = 0.05; t05 = tinv(1-alpha/2,n-p); int0 = [theta(1) - S*t05*sqrt(F(1,1)) theta(1) + S*t05*sqrt(F(1,1))] int1 = [theta(2) - S*t05*sqrt(F(2,2)) theta(2) + S*t05*sqrt(F(2,2))] int2 = [theta(3) - S*t05*sqrt(F(3,3)) theta(3) + S*t05*sqrt(F(3,3))] % % Figure 1: Data % Figure 2: Model fit % Figure 3: Residuals % figure(1) plot(x,y,'bo','linewidth',3) axis([1 9 144 158]) set(gca,'Xtick',[1 2 3 4 5 6 7 8 9]) set(gca,'XTickLabel',{'1','2','3','4','5','6','7','8','9'}) set(gca,'Fontsize',[22]); xlabel('Measurement Occasion') ylabel('Height (cm)') %print -depsc height_data figure(2) plot(x,y,'bo',x,Model,'k-','linewidth',3) axis([1 9 144 158]) set(gca,'Xtick',[1 2 3 4 5 6 7 8 9]) set(gca,'XTickLabel',{'1','2','3','4','5','6','7','8','9'}) set(gca,'Fontsize',[22]); legend('Data','Model','Location','SouthEast') xlabel('Measurement Occasion') ylabel('Height (cm)') %print -depsc height figure(3) plot(x,R,'bo',x,0*x,'k-','linewidth',3) axis([1 9 -0.8 0.8]) set(gca,'Fontsize',[22]); set(gca,'Xtick',[1 2 3 4 5 6 7 8 9]) set(gca,'XTickLabel',{'1','2','3','4','5','6','7','8','9'}) xlabel('Measurement Occasion') ylabel('Residual (cm)') %print -depsc height_residuals