Y=[1.7;.9;1.8];
X=[1,4.1+4+4.1;
    1,14.2;
    1,7+6.3];
hold off

plot3(Y(1),Y(2),Y(3),'*')
grid on
hold on 
plot3([0,X(1,1)],[0,X(2,1)],[0,X(3,1)])
plot3([0,X(1,2)],[0,X(2,2)],[0,X(3,2)])
xlabel('obs 1')
ylabel('obs 2')
zlabel('obs 3')
title('Marker height and width in the observation space and  the X vectors that will define the expectation plane')

%Plot the data in the observation space and the X vectors which will define the expectation plane


%%


% same plot but rescaling the intercept vector to make it more visible
Y=[1.7;.9;1.8];
X=[10,4.1+4+4.1;
    10,14.2;
    10,7+6.3];

figure(2)
plot3(Y(1),Y(2),Y(3),'*')
grid on
hold on 
plot3([0,X(1,1)],[0,X(2,1)],[0,X(3,1)])
plot3([0,X(1,2)],[0,X(2,2)],[0,X(3,2)])
xlabel('obs 1')
ylabel('obs 2')
zlabel('obs 3')
title('Marker height and width in the observation space and  the X vectors that will define the expectation plane')


%% Show The expectation plane


Y=[1.7;.9;1.8];
X=[1,4.1+4+4.1;
    1,14.2;
    1,7+6.3];


% evaluate the expectation surface over a grid of values for the two
% parameters
clear        zmat
beta1=-5:1:10;
beta2=-2:.5:2;
for lp=1:length(beta1)
    for kp=1:length(beta2)
        zmat(lp,kp,:)=X*[beta1(lp);beta2(kp)];
    end
end



% plot the expectation surface
hold off 
surf(zmat(:,:,1),zmat(:,:,2),zmat(:,:,3))

% rescale the plot
xlim([0,5])
ylim([0,5])
zlim([0,5])

% add the observation to the plot

hold on
plot3(Y(1),Y(2),Y(3),'*')
grid on
hold off


% now spin and  rotate the plot


%% QR
[Q,R]=qr(X)
N=3
p=2

Q1=Q(:,1:p);
Q2=Q(:,p+1:end);
R1=R(1:p,:);

w=Q'*Y;
w1=Q1'*Y;
hatY=Q1*w1;



% plot the expectation surface
hold off 
surf(zmat(:,:,1),zmat(:,:,2),zmat(:,:,3))


% rescale the plot
xlim([0,3])
ylim([0,3])
zlim([0,3])
% add the value of hatY
hold on 
plot3([Y(1),hatY(1)],[Y(2),hatY(2)],[Y(3),hatY(3)],'*-b')
plot3(hatY(1),hatY(2),hatY(3),'*r')

hold off

%%


% estimate the model parameters

hat_beta = R1\w1
regress(Y,X)


% residual error estimate:

Z=Y-hatY;
U=Q'*Z;
U1=U(1:p,:);
U2=U(p+1:end,:);

s2=U2'*U2/(N-p)
Z'*Z/(N-p)


figure
%%CI is
w=0:.1:2*pi+.05;
CI=repmat(hat_beta,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)]
plot(CI(1,:),CI(2,:),'r',hat_beta(1),hat_beta(2),'*')

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