wolffd@0
|
1 function [var, U, lambda] = ppca(x, ppca_dim)
|
wolffd@0
|
2 %PPCA Probabilistic Principal Components Analysis
|
wolffd@0
|
3 %
|
wolffd@0
|
4 % Description
|
wolffd@0
|
5 % [VAR, U, LAMBDA] = PPCA(X, PPCA_DIM) computes the principal
|
wolffd@0
|
6 % component subspace U of dimension PPCA_DIM using a centred covariance
|
wolffd@0
|
7 % matrix X. The variable VAR contains the off-subspace variance (which
|
wolffd@0
|
8 % is assumed to be spherical), while the vector LAMBDA contains the
|
wolffd@0
|
9 % variances of each of the principal components. This is computed
|
wolffd@0
|
10 % using the eigenvalue and eigenvector decomposition of X.
|
wolffd@0
|
11 %
|
wolffd@0
|
12 % See also
|
wolffd@0
|
13 % EIGDEC, PCA
|
wolffd@0
|
14 %
|
wolffd@0
|
15
|
wolffd@0
|
16 % Copyright (c) Ian T Nabney (1996-2001)
|
wolffd@0
|
17
|
wolffd@0
|
18
|
wolffd@0
|
19 if ppca_dim ~= round(ppca_dim) | ppca_dim < 1 | ppca_dim > size(x, 2)
|
wolffd@0
|
20 error('Number of PCs must be integer, >0, < dim');
|
wolffd@0
|
21 end
|
wolffd@0
|
22
|
wolffd@0
|
23 [ndata, data_dim] = size(x);
|
wolffd@0
|
24 % Assumes that x is centred and responsibility weighted
|
wolffd@0
|
25 % covariance matrix
|
wolffd@0
|
26 [l Utemp] = eigdec(x, data_dim);
|
wolffd@0
|
27 % Zero any negative eigenvalues (caused by rounding)
|
wolffd@0
|
28 l(l<0) = 0;
|
wolffd@0
|
29 % Now compute the sigma squared values for all possible values
|
wolffd@0
|
30 % of q
|
wolffd@0
|
31 s2_temp = cumsum(l(end:-1:1))./[1:data_dim]';
|
wolffd@0
|
32 % If necessary, reduce the value of q so that var is at least
|
wolffd@0
|
33 % eps * largest eigenvalue
|
wolffd@0
|
34 q_temp = min([ppca_dim; data_dim-min(find(s2_temp/l(1) > eps))]);
|
wolffd@0
|
35 if q_temp ~= ppca_dim
|
wolffd@0
|
36 wstringpart = 'Covariance matrix ill-conditioned: extracted';
|
wolffd@0
|
37 wstring = sprintf('%s %d/%d PCs', ...
|
wolffd@0
|
38 wstringpart, q_temp, ppca_dim);
|
wolffd@0
|
39 warning(wstring);
|
wolffd@0
|
40 end
|
wolffd@0
|
41 if q_temp == 0
|
wolffd@0
|
42 % All the latent dimensions have disappeared, so we are
|
wolffd@0
|
43 % just left with the noise model
|
wolffd@0
|
44 var = l(1)/data_dim;
|
wolffd@0
|
45 lambda = var*ones(1, ppca_dim);
|
wolffd@0
|
46 else
|
wolffd@0
|
47 var = mean(l(q_temp+1:end));
|
wolffd@0
|
48 end
|
wolffd@0
|
49 U = Utemp(:, 1:q_temp);
|
wolffd@0
|
50 lambda(1:q_temp) = l(1:q_temp);
|
wolffd@0
|
51
|
wolffd@0
|
52
|