function [yhat,e] = plot_col3(X,y,ngrid)

%   This file creates a 3-D plot.  The plot contains (a) the plane
%   generated by the column space of X that is orthogonal to N(X');
%   (b) the vector y (denoted by *); (c) the vector y-hat (denoted
%   by o); and the vector e = y - y-hat (denoted by * with a circle).
%   
%   Certain variables need to be input
%
%   X: 3 x p design matix with rank 2
%   y: 3 x 1 vector of observed responses
%   ngrid:  number of grid points in x and y directions.  Try ngrid =5
%   or ngrid = 100 to see what happens.
%
clf
o=zeros(3,1);
P=X*inv(X'*X)*X';
yhat=P*y;
P1=P(:,1:2);
P2=P(:,3);
I3=eye(3);
E12=I3(:,1:2);
e3=I3(:,3);
C=pinv(P2-e3)*(E12-P1);
e=y-yhat;
m1=min([X y yhat e]');
m2=max([X y yhat e]');
minx=min(0,m1(1))*1.25;
miny=min(0,m1(2))*1.25;
maxx=max(0,m2(1))*1.25;
maxy=max(0,m2(2))*1.25;
delx=(maxx-minx)/ngrid;
dely=(maxy-miny)/ngrid;
xx=[minx:delx:maxx]';
yy=[miny:dely:maxy]';
nx=length(xx);
ny=length(yy);
Z=zeros(nx,ny);
U=[];
for i=1:nx
  for j=1:ny
  Z(i,j)=C*[xx(i) yy(j)]';
    u=[xx(i) yy(j) Z(i,j)]';
    U=[U u];
  end;
end;
m1=min([U y yhat e]');
m2=max([U y yhat e]');
minx=min(0,m1(1));
miny=min(0,m1(2));
minz=min(0,m1(3));
maxx=max(0,m2(1));
maxy=max(0,m2(2));
maxz=max(0,m2(3));
mesh(xx,yy,Z')
axis([minx,maxx,miny,maxy,minz,maxz])
hold on
m1x=[maxx 0 0]';
m1y=[0 maxy 0]';
m1z=[0 0 maxz]';
m2x=[minx 0 0]';
m2y=[0 miny 0]';
m2z=[0 0 minz]';
V=[o';m1x';m2x';o';m1y';m2y';o';m1z';m2z'];
xx=V(:,1);
yy=V(:,2);
zz=V(:,3);
plot3(xx,yy,zz)
V=[o y yhat o e]';
xx=V(:,1);
yy=V(:,2);
zz=V(:,3);
plot3(xx,yy,zz)
plot3(y(1),y(2),y(3),'*')
plot3(yhat(1),yhat(2),yhat(3),'o')
e=y-yhat;
plot3(e(1),e(2),e(3),'*')
plot3(e(1),e(2),e(3),'o')
xlabel('Axis 1')
ylabel('Axis 2')
zlabel('Axis 3')