f1=figure;
numpoints=3;

B= mvnrnd([0 0], [1 .8; .8 1], numpoints)';
meanb=mean(B,2)
B=B-repmat(meanb,1,numpoints);
plot(B(1,:),B(2,:),'.')
axis equal;

hold on;

[Vp,Dp]=eig(B*B')
[V,D]=eigsort(Vp,Dp);
plot([0 V(1,1)],[0,V(2,1)])


%% plot the points in the usual x,y space
for i=1:numpoints
plot([0 B(1,i),B(1,i)],[0, 0, B(2,i)], 'r')
axis equal;

text(B(1,i), B(2,i), num2str([B(1,i),B(2,i)]))
end


f2=figure;
plot([0 V(1,1)] ,[0 V(2,1)])
axis equal;

%%% plot the points in the PC space
hold on
plot(0,0,'gx');
for i=1:numpoints
  c=V'*B(:,i);  %% gives PCA representation of ith data point
  midpt=c(1)*V(:,1);   %% This finds the point along the first principal
                       %component direction where the datapoint would be
                       %projected to
  plot([0 midpt(1),B(1,i)],[0, midpt(2), B(2,i)], 'r')  
  text(B(1,i), B(2,i), num2str([c(1),c(2)]))
end


f3=figure;
imagesc(V(:,1),[-10 10]);colormap(gray),colorbar;
imagesc(V(:,2),[-10 10]);colormap(gray),colorbar;


for i=-20:2:20
  imagesc(meanb+i*V(:,1),[-10 10]); colormap(gray),colorbar;
  M((i+22)/2)=getframe;
  pause;
end

movie(M)


C=V'*B

Breconstructed=V*C;

imagesc(B(:,1),[-1 1]); colormap(gray); colorbar
figure; imagesc(V(:,1)*C(1,:),[-1 1]); colormap(gray); colorbar

figure;imagesc(Breconstructed,[-1 1]); colormap(gray); colorbar