wolffd@0: function [var, U, lambda] = ppca(x, ppca_dim)
wolffd@0: %PPCA	Probabilistic Principal Components Analysis
wolffd@0: %
wolffd@0: %	Description
wolffd@0: %	 [VAR, U, LAMBDA] = PPCA(X, PPCA_DIM) computes the principal
wolffd@0: %	component subspace U of dimension PPCA_DIM using a centred covariance
wolffd@0: %	matrix X. The variable VAR contains the off-subspace variance (which
wolffd@0: %	is assumed to be spherical), while the vector LAMBDA contains the
wolffd@0: %	variances of each of the principal components.  This is computed
wolffd@0: %	using the eigenvalue and eigenvector  decomposition of X.
wolffd@0: %
wolffd@0: %	See also
wolffd@0: %	EIGDEC, PCA
wolffd@0: %
wolffd@0: 
wolffd@0: %	Copyright (c) Ian T Nabney (1996-2001)
wolffd@0: 
wolffd@0: 
wolffd@0: if ppca_dim ~= round(ppca_dim) | ppca_dim < 1 | ppca_dim > size(x, 2)
wolffd@0:    error('Number of PCs must be integer, >0, < dim');
wolffd@0: end
wolffd@0: 
wolffd@0: [ndata, data_dim] = size(x);
wolffd@0: % Assumes that x is centred and responsibility weighted
wolffd@0: % covariance matrix
wolffd@0: [l Utemp] = eigdec(x, data_dim);
wolffd@0: % Zero any negative eigenvalues (caused by rounding)
wolffd@0: l(l<0) = 0;
wolffd@0: % Now compute the sigma squared values for all possible values
wolffd@0: % of q
wolffd@0: s2_temp = cumsum(l(end:-1:1))./[1:data_dim]';
wolffd@0: % If necessary, reduce the value of q so that var is at least
wolffd@0: % eps * largest eigenvalue
wolffd@0: q_temp = min([ppca_dim; data_dim-min(find(s2_temp/l(1) > eps))]);
wolffd@0: if q_temp ~= ppca_dim
wolffd@0:   wstringpart = 'Covariance matrix ill-conditioned: extracted';
wolffd@0:   wstring = sprintf('%s %d/%d PCs', ...
wolffd@0:       wstringpart, q_temp, ppca_dim);
wolffd@0:   warning(wstring);
wolffd@0: end
wolffd@0: if q_temp == 0
wolffd@0:   % All the latent dimensions have disappeared, so we are
wolffd@0:   % just left with the noise model
wolffd@0:   var = l(1)/data_dim;
wolffd@0:   lambda = var*ones(1, ppca_dim);
wolffd@0: else
wolffd@0:   var = mean(l(q_temp+1:end));
wolffd@0: end  
wolffd@0: U = Utemp(:, 1:q_temp);
wolffd@0: lambda(1:q_temp) = l(1:q_temp);
wolffd@0: 
wolffd@0: