344 lines
11 KiB
Matlab
344 lines
11 KiB
Matlab
function [E, D] = pcamat(vectors, firstEig, lastEig, s_interactive, ...
|
|
s_verbose);
|
|
%PCAMAT - Calculates the pca for data
|
|
%
|
|
% [E, D] = pcamat(vectors, firstEig, lastEig, ...
|
|
% interactive, verbose);
|
|
%
|
|
% Calculates the PCA matrices for given data (row) vectors. Returns
|
|
% the eigenvector (E) and diagonal eigenvalue (D) matrices containing the
|
|
% selected subspaces. Dimensionality reduction is controlled with
|
|
% the parameters 'firstEig' and 'lastEig' - but it can also be done
|
|
% interactively by setting parameter 'interactive' to 'on' or 'gui'.
|
|
%
|
|
% ARGUMENTS
|
|
%
|
|
% vectors Data in row vectors.
|
|
% firstEig Index of the largest eigenvalue to keep.
|
|
% Default is 1.
|
|
% lastEig Index of the smallest eigenvalue to keep.
|
|
% Default is equal to dimension of vectors.
|
|
% interactive Specify eigenvalues to keep interactively. Note that if
|
|
% you set 'interactive' to 'on' or 'gui' then the values
|
|
% for 'firstEig' and 'lastEig' will be ignored, but they
|
|
% still have to be entered. If the value is 'gui' then the
|
|
% same graphical user interface as in FASTICAG will be
|
|
% used. Default is 'off'.
|
|
% verbose Default is 'on'.
|
|
%
|
|
%
|
|
% EXAMPLE
|
|
% [E, D] = pcamat(vectors);
|
|
%
|
|
% Note
|
|
% The eigenvalues and eigenvectors returned by PCAMAT are not sorted.
|
|
%
|
|
% This function is needed by FASTICA and FASTICAG
|
|
|
|
% For historical reasons this version does not sort the eigenvalues or
|
|
% the eigen vectors in any ways. Therefore neither does the FASTICA or
|
|
% FASTICAG. Generally it seams that the components returned from
|
|
% whitening is almost in reversed order. (That means, they usually are,
|
|
% but sometime they are not - depends on the EIG-command of matlab.)
|
|
|
|
% 16.6.2000
|
|
% Hugo Gävert
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% Default values:
|
|
if nargin < 5, s_verbose = 'on'; end
|
|
if nargin < 4, s_interactive = 'off'; end
|
|
if nargin < 3, lastEig = size(vectors, 1); end
|
|
if nargin < 2, firstEig = 1; end
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% Check the optional parameters;
|
|
if strcmp(lower(s_verbose), 'on')
|
|
b_verbose = 1;
|
|
elseif strcmp(lower(s_verbose), 'off')
|
|
b_verbose = 0;
|
|
else
|
|
error(sprintf('Illegal value [ %s ] for parameter: ''verbose''\n', s_verbose));
|
|
end
|
|
|
|
if strcmp(lower(s_interactive), 'on')
|
|
b_interactive = 1;
|
|
elseif strcmp(lower(s_interactive), 'off')
|
|
b_interactive = 0;
|
|
elseif strcmp(lower(s_interactive), 'gui')
|
|
b_interactive = 2;
|
|
else
|
|
error(sprintf('Illegal value [ %s ] for parameter: ''interactive''\n', ...
|
|
s_interactive));
|
|
end
|
|
|
|
oldDimension = size (vectors, 1);
|
|
if ~(b_interactive)
|
|
if lastEig < 1 | lastEig > oldDimension
|
|
error(sprintf('Illegal value [ %d ] for parameter: ''lastEig''\n', lastEig));
|
|
end
|
|
if firstEig < 1 | firstEig > lastEig
|
|
error(sprintf('Illegal value [ %d ] for parameter: ''firstEig''\n', firstEig));
|
|
end
|
|
end
|
|
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% Calculate PCA
|
|
|
|
% Calculate the covariance matrix.
|
|
if b_verbose, fprintf ('Calculating covariance...\n'); end
|
|
covarianceMatrix = cov(vectors')*(size(vectors',1)-1)/size(vectors',1);
|
|
|
|
maxLastEig = rank(covarianceMatrix, 1e-9);
|
|
|
|
% Calculate the eigenvalues and eigenvectors of covariance matrix.
|
|
[E, D] = eig(covarianceMatrix);
|
|
|
|
% Sort the eigenvalues - decending.
|
|
eigenvalues = flipud(sort(diag(D)));
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% Interactive part - command-line
|
|
if b_interactive == 1
|
|
|
|
% Show the eigenvalues to the user
|
|
hndl_win=figure;
|
|
bar(eigenvalues);
|
|
title('Eigenvalues');
|
|
|
|
% ask the range from the user...
|
|
% ... and keep on asking until the range is valid :-)
|
|
areValuesOK=0;
|
|
while areValuesOK == 0
|
|
firstEig = input('The index of the largest eigenvalue to keep? (1) ');
|
|
lastEig = input(['The index of the smallest eigenvalue to keep? (' ...
|
|
int2str(oldDimension) ') ']);
|
|
% Check the new values...
|
|
% if they are empty then use default values
|
|
if isempty(firstEig), firstEig = 1;end
|
|
if isempty(lastEig), lastEig = oldDimension;end
|
|
% Check that the entered values are within the range
|
|
areValuesOK = 1;
|
|
if lastEig < 1 | lastEig > oldDimension
|
|
fprintf('Illegal number for the last eigenvalue.\n');
|
|
areValuesOK = 0;
|
|
end
|
|
if firstEig < 1 | firstEig > lastEig
|
|
fprintf('Illegal number for the first eigenvalue.\n');
|
|
areValuesOK = 0;
|
|
end
|
|
end
|
|
% close the window
|
|
close(hndl_win);
|
|
end
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% Interactive part - GUI
|
|
if b_interactive == 2
|
|
|
|
% Show the eigenvalues to the user
|
|
hndl_win = figure('Color',[0.8 0.8 0.8], ...
|
|
'PaperType','a4letter', ...
|
|
'Units', 'normalized', ...
|
|
'Name', 'FastICA: Reduce dimension', ...
|
|
'NumberTitle','off', ...
|
|
'Tag', 'f_eig');
|
|
h_frame = uicontrol('Parent', hndl_win, ...
|
|
'BackgroundColor',[0.701961 0.701961 0.701961], ...
|
|
'Units', 'normalized', ...
|
|
'Position',[0.13 0.05 0.775 0.17], ...
|
|
'Style','frame', ...
|
|
'Tag','f_frame');
|
|
|
|
b = uicontrol('Parent',hndl_win, ...
|
|
'Units','normalized', ...
|
|
'BackgroundColor',[0.701961 0.701961 0.701961], ...
|
|
'HorizontalAlignment','left', ...
|
|
'Position',[0.142415 0.0949436 0.712077 0.108507], ...
|
|
'String','Give the indices of the largest and smallest eigenvalues of the covariance matrix to be included in the reduced data.', ...
|
|
'Style','text', ...
|
|
'Tag','StaticText1');
|
|
e_first = uicontrol('Parent',hndl_win, ...
|
|
'Units','normalized', ...
|
|
'Callback',[ ...
|
|
'f=round(str2num(get(gcbo, ''String'')));' ...
|
|
'if (f < 1), f=1; end;' ...
|
|
'l=str2num(get(findobj(''Tag'',''e_last''), ''String''));' ...
|
|
'if (f > l), f=l; end;' ...
|
|
'set(gcbo, ''String'', int2str(f));' ...
|
|
], ...
|
|
'BackgroundColor',[1 1 1], ...
|
|
'HorizontalAlignment','right', ...
|
|
'Position',[0.284831 0.0678168 0.12207 0.0542535], ...
|
|
'Style','edit', ...
|
|
'String', '1', ...
|
|
'Tag','e_first');
|
|
b = uicontrol('Parent',hndl_win, ...
|
|
'Units','normalized', ...
|
|
'BackgroundColor',[0.701961 0.701961 0.701961], ...
|
|
'HorizontalAlignment','left', ...
|
|
'Position',[0.142415 0.0678168 0.12207 0.0542535], ...
|
|
'String','Range from', ...
|
|
'Style','text', ...
|
|
'Tag','StaticText2');
|
|
e_last = uicontrol('Parent',hndl_win, ...
|
|
'Units','normalized', ...
|
|
'Callback',[ ...
|
|
'l=round(str2num(get(gcbo, ''String'')));' ...
|
|
'lmax = get(gcbo, ''UserData'');' ...
|
|
'if (l > lmax), l=lmax; fprintf([''The selected value was too large, or the selected eigenvalues were close to zero\n'']); end;' ...
|
|
'f=str2num(get(findobj(''Tag'',''e_first''), ''String''));' ...
|
|
'if (l < f), l=f; end;' ...
|
|
'set(gcbo, ''String'', int2str(l));' ...
|
|
], ...
|
|
'BackgroundColor',[1 1 1], ...
|
|
'HorizontalAlignment','right', ...
|
|
'Position',[0.467936 0.0678168 0.12207 0.0542535], ...
|
|
'Style','edit', ...
|
|
'String', int2str(maxLastEig), ...
|
|
'UserData', maxLastEig, ...
|
|
'Tag','e_last');
|
|
% in the first version oldDimension was used instead of
|
|
% maxLastEig, but since the program would automatically
|
|
% drop the eigenvalues afte maxLastEig...
|
|
b = uicontrol('Parent',hndl_win, ...
|
|
'Units','normalized', ...
|
|
'BackgroundColor',[0.701961 0.701961 0.701961], ...
|
|
'HorizontalAlignment','left', ...
|
|
'Position',[0.427246 0.0678168 0.0406901 0.0542535], ...
|
|
'String','to', ...
|
|
'Style','text', ...
|
|
'Tag','StaticText3');
|
|
b = uicontrol('Parent',hndl_win, ...
|
|
'Units','normalized', ...
|
|
'Callback','uiresume(gcbf)', ...
|
|
'Position',[0.630697 0.0678168 0.12207 0.0542535], ...
|
|
'String','OK', ...
|
|
'Tag','Pushbutton1');
|
|
b = uicontrol('Parent',hndl_win, ...
|
|
'Units','normalized', ...
|
|
'Callback',[ ...
|
|
'gui_help(''pcamat'');' ...
|
|
], ...
|
|
'Position',[0.767008 0.0678168 0.12207 0.0542535], ...
|
|
'String','Help', ...
|
|
'Tag','Pushbutton2');
|
|
|
|
h_axes = axes('Position' ,[0.13 0.3 0.775 0.6]);
|
|
set(hndl_win, 'currentaxes',h_axes);
|
|
bar(eigenvalues);
|
|
title('Eigenvalues');
|
|
|
|
uiwait(hndl_win);
|
|
firstEig = str2num(get(e_first, 'String'));
|
|
lastEig = str2num(get(e_last, 'String'));
|
|
|
|
% close the window
|
|
close(hndl_win);
|
|
end
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% See if the user has reduced the dimension enought
|
|
|
|
if lastEig > maxLastEig
|
|
lastEig = maxLastEig;
|
|
if b_verbose
|
|
fprintf('Dimension reduced to %d due to the singularity of covariance matrix\n',...
|
|
lastEig-firstEig+1);
|
|
end
|
|
else
|
|
% Reduce the dimensionality of the problem.
|
|
if b_verbose
|
|
if oldDimension == (lastEig - firstEig + 1)
|
|
fprintf ('Dimension not reduced.\n');
|
|
else
|
|
fprintf ('Reducing dimension...\n');
|
|
end
|
|
end
|
|
end
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% Drop the smaller eigenvalues
|
|
if lastEig < oldDimension
|
|
lowerLimitValue = (eigenvalues(lastEig) + eigenvalues(lastEig + 1)) / 2;
|
|
else
|
|
lowerLimitValue = eigenvalues(oldDimension) - 1;
|
|
end
|
|
|
|
lowerColumns = diag(D) > lowerLimitValue;
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% Drop the larger eigenvalues
|
|
if firstEig > 1
|
|
higherLimitValue = (eigenvalues(firstEig - 1) + eigenvalues(firstEig)) / 2;
|
|
else
|
|
higherLimitValue = eigenvalues(1) + 1;
|
|
end
|
|
higherColumns = diag(D) < higherLimitValue;
|
|
|
|
% Combine the results from above
|
|
selectedColumns = lowerColumns & higherColumns;
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% print some info for the user
|
|
if b_verbose
|
|
fprintf ('Selected [ %d ] dimensions.\n', sum (selectedColumns));
|
|
end
|
|
if sum (selectedColumns) ~= (lastEig - firstEig + 1),
|
|
error ('Selected a wrong number of dimensions.');
|
|
end
|
|
|
|
if b_verbose
|
|
fprintf ('Smallest remaining (non-zero) eigenvalue [ %g ]\n', eigenvalues(lastEig));
|
|
fprintf ('Largest remaining (non-zero) eigenvalue [ %g ]\n', eigenvalues(firstEig));
|
|
fprintf ('Sum of removed eigenvalues [ %g ]\n', sum(diag(D) .* ...
|
|
(~selectedColumns)));
|
|
end
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% Select the colums which correspond to the desired range
|
|
% of eigenvalues.
|
|
E = selcol (E, selectedColumns);
|
|
D = selcol (selcol (D, selectedColumns)', selectedColumns);
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% Some more information
|
|
if b_verbose
|
|
sumAll=sum(eigenvalues);
|
|
sumUsed=sum(diag(D));
|
|
retained = (sumUsed / sumAll) * 100;
|
|
fprintf('[ %g ] %% of (non-zero) eigenvalues retained.\n', retained);
|
|
end
|
|
endfunction
|
|
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
|
|
% Moved to "selcol.m"
|
|
|
|
%function newMatrix = selcol(oldMatrix, maskVector);
|
|
|
|
% newMatrix = selcol(oldMatrix, maskVector);
|
|
|
|
% Selects the columns of the matrix that marked by one in the given vector.
|
|
% The maskVector is a column vector.
|
|
|
|
% 15.3.1998
|
|
|
|
%if size(maskVector, 1) ~= size(oldMatrix, 2),
|
|
% error ('The mask vector and matrix are of uncompatible size.');
|
|
%end
|
|
|
|
%numTaken = 0;
|
|
|
|
%for i = 1 : size (maskVector, 1),
|
|
% if maskVector(i, 1) == 1,
|
|
% takingMask(1, numTaken + 1) = i;
|
|
% numTaken = numTaken + 1;
|
|
% end
|
|
%end
|
|
|
|
%newMatrix = oldMatrix(:, takingMask);
|
|
%endfunction
|