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