function [mu,V,X,c,seed] = mnist_pca(n,nums,seed);
% function [mu,V,X] = mnist_pca(n,nums,probobs);
% Returns PCA reduced MNIST, with mu+V*X as approximation to images (one per col)
% n = number of dimensions in PCA.
% nums = number of images of each digit to use (between 1 and 100)

if nargin < 1 || isempty(n),       n       = 50; end
if nargin < 2 || isempty(nums),    nums    = 50; end
if nargin < 3 || isempty(probobs), probobs = 1; end
if nargin < 4 || isempty(seed),    
  seed = rand('twister'); 
  seed2 = [];
else
  seed2 = rand('twister');
  rand('twister',seed);
end

load mnist
data = D;
m = 10*nums;
c = reshape(repmat(1:10,nums,1),[1 m]); % class labels
D = zeros(784,m);
for k = 1:10
  rr = randperm(100);
  D(:,(k-1)*nums+(1:nums)) = reshape(permute(reshape(data(rr(1:nums),:,k),...
        [nums 28 28]),[3 2 1]),[784 nums])/255;
end
mu = mean(D,2);
[V S Y] = svd(D-mu(:,ones(1,m)),0);
V = V*S(:,1:n)/sqrt(m);
X = (Y(:,1:n)'*sqrt(m))'; % data

if ~isempty(seed2),
  rand('twister',seed2);
end


