clc
clear
% Bates and Watts table 1.16 lipoprotein tracer concentration versus time

time=[.5   ,1   ,1.5,2    ,3   ,4   ,5   ,6   ,7   ,8   ,9 ,10];
data=[46.1,25.9,17  ,12.10,7.22,4.51,3.19,2.4,1.82,1.41 ,1 ,.94];

plot(time,data,'.')


% we will use a subset of the data

X=[1,3, 7]'
Y=[25.9,7.22,1.82]'


plot3(Y(1),Y(2),Y(3),'*')


cd('/Volumes/iamdavecampbell/stat890-2009')


%the model for X
clear zmat

theta=-.5:.25:10;
for lp=1:length(theta)
    zmat(lp,:)=lipomodel(X,theta(lp));
end





% plot the expectation surface

plot3(zmat(:,1),zmat(:,2),zmat(:,3),'.-')



% add the observation to the plot

hold on
plot3(Y(1),Y(2),Y(3),'*k','MarkerSize',20,'LineWidth',3)
grid on
hold off

xlabel('\fontsize{18}obs #1')
ylabel('\fontsize{18}obs #2')
zlabel('\fontsize{18}obs #3')
title('\fontsize{20}lipoprotein observation and the expectation surface in the observation space')
% dots in expectation surface are evenly spaced at intervals of .25 in
% theta

% rescale the plot to zoom in on the features
xlim([0,100])
ylim([0,100])
zlim([0,100])


xlim([0,40])
ylim([0,10])
zlim([0,3])



%%
% 2 dimensional case

cd('/Volumes/iamdavecampbell/stat890-2009')


%the model for X
clear     zmat
theta2=0:.25:10;
theta1=10:10:100;
for lp=1:length(theta1)
    for kp=1:length(theta2)
        zmat(lp,kp,:)=lipomodel(X,[theta1(lp),theta2(kp)]);
    end
end


% plot the expectation surface
figure(1)
plot3(zmat(:,:,1)',zmat(:,:,2)',zmat(:,:,3)','.-')
%fixed by adding transpost to this plot3


% add the observation to the plot

hold on
plot3(Y(1),Y(2),Y(3),'*k','MarkerSize',20,'LineWidth',3)
grid on
hold off

xlabel('\fontsize{18}obs #1')
ylabel('\fontsize{18}obs #2')
zlabel('\fontsize{18}obs #3')
title('\fontsize{20}lipoprotein observation and the expectation surface in the observation space')
% dots in expectation surface are evenly spaced at intervals of .25 in
% theta

% rescale the plot to zoom in on the features
xlim([0,40])
ylim([0,10])
zlim([0,3])


%%  Gauss-Newton Estimation
%starting values here were chosen so that we assume the initial value is
%X(0)=100 and the rate parameter is somewhat arbitrarily chosen
theta=[100,1.5]';%
% we need the derivative function  lipoproteinmodeldfdp
p=2
N=length(Y)

niter=100;
lp=1


tol = 10^-4
delta=100
while lp<niter && abs(max(delta))>tol
    etahat(:,lp) = lipomodel(X,theta(:,lp));
    V= lipomodeldfdp(X,theta(:,lp));
    [Q, R]=qr(V);
    Q1=Q(:,1:p);
    R1=R(1:p,:);
    z = Y-etahat(:,lp);
    w1=Q1'*z;
    delta=R1\w1;
    theta(:,lp+1)=theta(:,lp)+delta;
    lp=lp+1;
end


% Plot the resultsand the expectation surface


%the model for X
clear     zmat
theta1=0:1:100;
theta2=0:.1:10;
for lp=1:length(theta1)
    for kp=1:length(theta2)
        zmat(lp,kp,:)=lipomodel(X,[theta1(lp),theta2(kp)]);
    end
end



figure(1)
surf(zmat(:,:,1),zmat(:,:,2),zmat(:,:,3))

hold on
plot3(Y(1),Y(2),Y(3),'*k','MarkerSize',20,'LineWidth',3)
grid on


plot3(etahat(1,:),etahat(2,:),etahat(3,:),'*-m','LineWidth',3)
plot3(etahat(1,end),etahat(2,end),etahat(3,end),'*r','MarkerSize',20,'LineWidth',3)


hold off

% note that the first step went way too far in the specified direction


xlabel('\fontsize{18}obs #1')
ylabel('\fontsize{18}obs #2')
zlabel('\fontsize{18}obs #3')
title('\fontsize{20}lipoprotein observations and response surface estimate')











% rescale the plot to zoom in on the features
xlim([0,40])
ylim([0,10])
zlim([0,5])



%%  Gauss Newton Estimation with reduced step size

%starting values here were chosen so that we assume the initial value is
%X(0)=100 and the rate parameter is somewhat arbitrarily chosen
theta=[100,1.5]';%
% we need the derivative function  lipoproteinmodeldfdp


niter=100;
lp=1

tol = 10^-4
delta=100
while lp<niter && abs(max(delta))>tol
    % start with everything as we did it before
    
    etahat(:,lp) = lipomodel(X,theta(:,lp));
    V= lipomodeldfdp(X,theta(:,lp));
    [Q, R]=qr(V);
    Q1=Q(:,1:p);
    R1=R(1:p,:);
    z = Y-etahat(:,lp);
    w1=Q1'*z;
    delta=R1\w1;
    discount = 1;
    theta(:,lp+1)=theta(:,lp)+delta*discount;
    
    % now see if the step is too far
    etahat(:,lp+1) = lipomodel(X,theta(:,lp+1));
    V= lipomodeldfdp(X,theta(:,lp+1));
    [Q, R]=qr(V);
    Q1=Q(:,1:p);
    R1=R(1:p,:);
    ztemp = Y-etahat(:,lp+1);
          
    while z'*z-ztemp'*ztemp<0   % while the step is too large
        delta=delta/2;
        theta(:,lp+1)=theta(:,lp)+delta*discount;
        etahat(:,lp+1) = lipomodel(X,theta(:,lp+1));
        V= lipomodeldfdp(X,theta(:,lp+1));
        [Q, R]=qr(V);
        Q1=Q(:,1:p);
        R1=R(1:p,:);
        ztemp = Y-etahat(:,lp+1);
        z=ztemp; % I will need the final residual value later so I want to keep it here
    % the last values computed in here will be the ones that we want to keep
    end
    lp=lp+1;
end





%the model for X
clear     zmat
theta1=0:1:100;
theta2=0:.1:10;
for lp=1:length(theta1)
    for kp=1:length(theta2)
        zmat(lp,kp,:)=lipomodel(X,[theta1(lp),theta2(kp)]);
    end
end



figure(1)
surf(zmat(:,:,1),zmat(:,:,2),zmat(:,:,3))

hold on
plot3(Y(1),Y(2),Y(3),'*k','MarkerSize',20,'LineWidth',3)
grid on


plot3(etahat(1,:),etahat(2,:),etahat(3,:),'-r','LineWidth',3)
plot3(etahat(1,end),etahat(2,end),etahat(3,end),'*r','MarkerSize',20,'LineWidth',3)


hold off

% note that the first step went way too far in the specified direction


xlabel('\fontsize{18}obs #1')
ylabel('\fontsize{18}obs #2')
zlabel('\fontsize{18}obs #3')
title('\fontsize{20}lipoprotein observations and response surface estimate with adjusted Gauss Newton Steps')











% rescale the plot to zoom in on the features
xlim([0,40])
ylim([0,10])
zlim([0,5])




%% make a joint confidence interval for theta1 and theta2

hattheta=theta(:,end)  % all we want is the final parameter estimate.  I only kept the intermediate values to show how it moves

w=0:.1:2*pi+.05;
s2=z'*z/(N-p)
CI=repmat(hattheta,1,length(w))+sqrt(p*s2*finv(.95,p,N-p))*R1\[cos(w);sin(w)]


% note that R1\[cos(w);sin(w)]  == R1^-1*[cos(w);sin(w)] == inv(R1)*[cos(w);sin(w)]
figure(11)
plot(CI(1,:),CI(2,:),'r',hattheta(1),hattheta(2),'*r')
xlabel('\fontsize{18}\theta_1')
ylabel('\fontsize{18}\theta_1')
title('\fontsize{20}CI for lipoprotein 3 observation model')

% label the plot
xlabel('\beta_1')
ylabel('\beta_2')



%% plot the confidence interval for the expected response
clear etaCI
x=0:.25:10;
s2=z'*z/(N-p);
for lp=1:length(x)
    V= lipomodeldfdp(x(lp),hattheta);
    SE=sqrt(s2*p*finv(.975,p,N-p)*(  V/R1 )*(  V/R1 )');
    etaCI(2,lp)=lipomodel(x(lp),hattheta);
    etaCI(1,lp)=etaCI(2,lp)- SE;
    etaCI(3,lp)=etaCI(2,lp)+ SE;
    
end


%the model for X
clear     zmat
theta1=0:1:100;
theta2=0:.1:10;
for lp=1:length(theta1)
    for kp=1:length(theta2)
        zmat(lp,kp,:)=lipomodel(X,[theta1(lp),theta2(kp)]);
    end
end



figure(12)


plot(x,etaCI,'r')
hold on
plot(X,Y,'*r')
title('\fontsize{20}3 Observations and the interval estimate for the function')





%%










































clear etahat 
X=time'%[.5   ,1   ,1.5,2    ,3   ,4   ,5   ,6   ,7   ,8   ,9 ,10];
Y=data'%=[46.1,25.9,17  ,12.10,7.22,4.51,3.19,2.4,1.82,1.41 ,1 ,.94];
%  Gauss Newton Estimation with reduced step size

%starting values here were chosen so that we assume the initial value is
%X(0)=100 and the rate parameter is somewhat arbitrarily chosen
theta=[100,1.5]';%
% we need the derivative function  lipoproteinmodeldfdp

N=length(Y)
p=2
niter=100;
lp=1

tol = 10^-4
delta=100
while lp<niter && abs(max(delta))>tol
    % start with everything as we did it before
    
    etahat(:,lp) = lipomodel(X,theta(:,lp));
    V= lipomodeldfdp(X,theta(:,lp));
    [Q, R]=qr(V);
    Q1=Q(:,1:p);
    R1=R(1:p,:);
    z = Y-etahat(:,lp);
    w1=Q1'*z;
    delta=R1\w1;
    discount = 1;
    theta(:,lp+1)=theta(:,lp)+delta*discount;
    
    % now see if the step is too far
    etahat(:,lp+1) = lipomodel(X,theta(:,lp+1));
    V= lipomodeldfdp(X,theta(:,lp+1));
    [Q, R]=qr(V);
    Q1=Q(:,1:p);
    R1=R(1:p,:);
    ztemp = Y-etahat(:,lp+1);
          
    while z'*z-ztemp'*ztemp<0   % while the step is too large
        delta=delta/2;
        theta(:,lp+1)=theta(:,lp)+delta*discount;
        etahat(:,lp+1) = lipomodel(X,theta(:,lp+1));
        V= lipomodeldfdp(X,theta(:,lp+1));
        [Q, R]=qr(V);
        Q1=Q(:,1:p);
        R1=R(1:p,:);
        ztemp = Y-etahat(:,lp+1);
        z=ztemp; % I will need the final residual value later so I want to keep it here
    % the last values computed in here will be the ones that we want to keep
    end
    lp=lp+1;
end



% make a joint confidence interval for theta1 and theta2

hattheta=theta(:,end)  % all we want is the final parameter estimate.  I only kept the intermediate values to show how it moves

w=0:.1:2*pi+.05;
s2=z'*z/(N-p);
CI=repmat(hattheta,1,length(w))+sqrt(p*s2*finv(.95,p,N-p))*R1\[cos(w);sin(w)]


% note that R1\[cos(w);sin(w)]  == R1^-1*[cos(w);sin(w)] == inv(R1)*[cos(w);sin(w)]
figure(11)
hold on
plot(CI(1,:),CI(2,:),'k',hattheta(1),hattheta(2),'*k')
xlabel('\fontsize{18}\theta_1')
ylabel('\fontsize{18}\theta_1')
title('\fontsize{20}CI for lipoprotein 3 observation model')

% label the plot
xlabel('\beta_1')
ylabel('\beta_2')

hold off

%% plot the confidence interval for the expected response
clear etaCI
x=0:.25:10;
s2=z'*z/(N-p);
for lp=1:length(x)
    V= lipomodeldfdp(x(lp),hattheta);
    SE=sqrt(s2*p*finv(.975,p,N-p)*(  V/R1 )*(  V/R1 )');
    etaCI(2,lp)=lipomodel(x(lp),hattheta);
    etaCI(1,lp)=etaCI(2,lp)- SE;
    etaCI(3,lp)=etaCI(2,lp)+ SE;
    
end





figure(12)

hold on
plot(x,etaCI,'k')
hold on
plot(X,Y,'*k')
title(strcat('\fontsize{20}',num2str(N),' Observations and the interval estimate for the function'))


hold off


%%  Using Matlab NLS function
theta=[100,1.5]'

[hattheta_NLS,z_NLS,J,COVB,mse] = nlinfit(X,Y,@lipomodel_Xlast,theta);
[z'*z/(N-p),mse] % these are the same!
[hattheta,hattheta_NLS]% values are the same using straight QR Gauss-Newton or the built in NLS function (which also uses Gauss Newton and QR.

CI_NLS =  nlparci(hattheta_NLS,z_NLS,'jacobian',J)

%% Get CIs for curves using

nlintool(X,Y,@lipomodel_Xlast,theta)


