function [var, U, lambda] = ppca(x, ppca_dim)
%PPCA	Probabilistic Principal Components Analysis%%	Description%	 [VAR, U, LAMBDA] = PPCA(X, PPCA_DIM) computes the principal%	component subspace U of dimension PPCA_DIM using a centred covariance%	matrix X. The variable VAR contains the off-subspace variance (which%	is assumed to be spherical), while the vector LAMBDA contains the%	variances of each of the principal components.  This is computed%	using the eigenvalue and eigenvector  decomposition of X.%%	See also%	EIGDEC, PCA%%	Copyright (c) Ian T Nabney (1996-2001)
if ppca_dim ~= round(ppca_dim) | ppca_dim < 1 | ppca_dim > size(x, 2)
   error('Number of PCs must be integer, >0, < dim');
end

[ndata, data_dim] = size(x);
% Assumes that x is centred and responsibility weighted
% covariance matrix
[l Utemp] = eigdec(x, data_dim);
% Zero any negative eigenvalues (caused by rounding)
l(l<0) = 0;
% Now compute the sigma squared values for all possible values
% of q
s2_temp = cumsum(l(end:-1:1))./[1:data_dim]';
% If necessary, reduce the value of q so that var is at least
% eps * largest eigenvalue
q_temp = min([ppca_dim; data_dim-min(find(s2_temp/l(1) > eps))]);
if q_temp ~= ppca_dim
  wstringpart = 'Covariance matrix ill-conditioned: extracted';
  wstring = sprintf('%s %d/%d PCs', ...
      wstringpart, q_temp, ppca_dim);
  warning(wstring);
end
if q_temp == 0
  % All the latent dimensions have disappeared, so we are
  % just left with the noise model
  var = l(1)/data_dim;
  lambda = var*ones(1, ppca_dim);
else
  var = mean(l(q_temp+1:end));
end  
U = Utemp(:, 1:q_temp);
lambda(1:q_temp) = l(1:q_temp);