Skip to content

GitLab

Commit 3fb9b89e authored by Saad Jbabdi's avatar Saad Jbabdi
Browse files

updates to ICA prac

parent 28de6ae6
Expand all Hide whitespace changes
Inline Side-by-side
independent_component_analysis/ica_prac.m
View file @ 3fb9b89e
......@@ -19,7 +19,7 @@
% we can optimise it) is what FastICA does
% %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1. Generate a points cloud in 2d
%% 1. Generate a points cloud in 2d
rng('default') % for replication
n = 400; % number of samples
d = 2; % dimensionality
......@@ -63,11 +63,11 @@ axis equal
% Centering (mean(data)=0 along each dimension)
C = eye(n) - ones(n,1)*ones(1,n)/n;
x_c = x*C;
% % optionally, we can use call 'mean'
% % A more efficient alternative is to simply call 'mean'
% x_c = x - mean(x,2);
% Pre-whitening (var(data)=1 along each dimension)
% This means that PCA becomes irrelevant: all direction have the same
% Pre-whitening (such that var(data)=1 along each dimension)
% This means that PCA becomes irrelevant: all directions have the same
% variance
[V,D] = eig(cov(x'));
Q = V*diag(1./sqrt(diag(D)))*V';
......@@ -95,12 +95,17 @@ axis equal
mean(x_c,2)
%%
%% quick check that the variance is 1 along all directions
addpath ./utils
G = @(x)(x.^2);
plot_cf(x_w,G)
figure
subplot(1,2,1), plot_cf(x_c,G)
legend({'data','variance'}), title('Before whitening');
subplot(1,2,2), plot_cf(x_w,G)
legend({'data','variance'}), title('After whitening');
%% random projection
%% random projections
figure
% 10 random projections
% (press any key to see the next projection)
......@@ -117,39 +122,20 @@ for i = 1:10
end
%% ???
t = linspace(0,10*pi,400);
s1 = sin(t) + .1*randn(size(t));
s2 = (mod(4*t,pi)) + .1*randn(size(t));
A = [1 .1;-0.3 1];
X = A*[s1;s2];
[X_w,~] = whiten(X);
W = plot_cost(X_w);
quiver(0,0,3*W(1,1),3*W(1,2),'linewidth',3,'color','k')
quiver(0,0,3*W(2,1),3*W(2,2),'linewidth',3,'color','k')
figure,
subplot(1,2,1)
plot(X_w')
subplot(1,2,2)
plot(((W*X_w))')
%% quick PCA
% You should be able to see that PCA picks up directions of high variance
% which do not in this case coincide with the interesting directions
% in this data
[V,~] = eig(cov(x'));
figure
% plot data
plot(x(1,:),x(2,:),'o'); hold on
% plot 1st PC direction (pos and neg)
% plot 1st PC direction in RED (pos and neg)
quiver(mean(x(1,:)),mean(x(2,:)),V(1,2),V(2,2),'linewidth',3,'color','r')
quiver(mean(x(1,:)),mean(x(2,:)),-V(1,2),-V(2,2),'linewidth',3,'color','r')
% plot 2nd PC direction (pos and neg)
% plot 2nd PC direction in GREEN (pos and neg)
quiver(mean(x(1,:)),mean(x(2,:)),V(1,1),V(2,1),'linewidth',3,'color','g')
quiver(mean(x(1,:)),mean(x(2,:)),-V(1,1),-V(2,1),'linewidth',3,'color','g')
axis equal
......@@ -162,7 +148,7 @@ grid on
% these are functions which, when applied to the data, 'highlight' the
% non-Gaussianities
% they also are supposed to approximate neg-entropy
G0 = @(x)(x.^2); % quadratic (do not use!)
G0 = @(x)(x.^2); % quadratic (do not use this!)
G1 = @(x)(1-exp(-x.^2)); % Gaussian --> integral(x*exp(-x^2/2)
G2 = @(x)(log(cosh(x))); % --> integral(tanh)
G3 = @(x)(x.^3); % kurtosis
......@@ -222,27 +208,33 @@ legend({'data','non-gaussianity','variance'})
%% Plot in de-whitened space
% generate bunch of points
[xx,yy] = pol2cart(angles,I/norm(I)*40);
ww=Q\[xx;yy];
% un-whiten with the whitening matrix Q
ww = Q\[xx;yy];
figure,hold on
plot(x_c(1,:),x_c(2,:),'o')
plot(ww(1,:),ww(2,:),'linewidth',3,'color','r')
quiver(0,0,V(1,2),V(2,2),'linewidth',3,'color','k')
quiver(0,0,-V(1,2),-V(2,2),'linewidth',3,'color','k')
axis equal
grid on
%% Compare with the FastICA toolbox :)
% Can be downloaded from https://research.ics.aalto.fi/ica/fastica/
addpath ./utils/FastICA_25
%% Compare to FastICA toolbox :)
% check the help to see how to change things like the contrast function g
% note that we are using x (fastica does pre-whitening internally)
[S, A, W] = fastica(x,'approach', 'symm','g','skew');
addpath ~/matlab/FastICA_25/
[W,~] = fastica(x,'approach', 'symm','g','skew');
W = W';
% I am not sure why I need to do the below
W = A';
% plot the results
figure
plot(x_c(1,:),x_c(2,:),'o'); hold on
% plot(w(1,:),w(2,:),'linewidth',3,'color','r')
quiver(0,0,W(1,1),W(1,2),'linewidth',3,'color','k')
quiver(0,0,-W(1,1),-W(1,2),'linewidth',3,'color','k')
quiver(0,0,W(2,1),W(2,2),'linewidth',3,'color','c')
......@@ -252,7 +244,7 @@ axis equal
%% Quick implementation of the FastICA algorithm
d = size(x,1);
W = orth(randn(d,d)) % Random initial weights
W = orth(randn(d,d)); % Random initial weights
delta = inf; % initialize delta (w_{i+1} - w_{i})
k = 0; % counting iteration
% maxiter = 1e4; % stop iteration if k = maxiter
......@@ -271,9 +263,9 @@ while delta > tolerance
Gp = 6 * Sk.^2; % derivative of G
W = (G * x_w') / n - bsxfun(@times,mean(Gp,2),W);
W = normRows(W);
W = W ./ sqrt(sum(W.^2,2));
% Decorrelate weights
% Decorrelate W's
[U, S, ~] = svd(W,'econ');
W = U * diag(1 ./ diag(S)) * U' * W;
......
independent_component_analysis/utils/FastICA_25/CVS/Entries 0 → 100644
View file @ 3fb9b89e
/Contents.m/1.7/Wed Oct 19 13:05:33 2005//
/demosig.m/1.2/Sat Apr 5 14:23:57 2003//
/dispsig.m/1.2/Sat Apr 5 14:23:57 2003//
/fastica.m/1.14/Wed Oct 19 13:05:34 2005//
/fasticag.m/1.5/Wed Oct 19 13:05:34 2005//
/fpica.m/1.7/Thu Jun 16 12:52:55 2005//
/gui_adv.m/1.4/Tue Jul 27 13:09:26 2004//
/gui_advc.m/1.3/Mon Sep 8 11:28:58 2003//
/gui_cb.m/1.5/Wed Sep 10 10:33:41 2003//
/gui_cg.m/1.2/Sat Apr 5 14:23:57 2003//
/gui_help.m/1.6/Wed Oct 19 13:05:34 2005//
/gui_l.m/1.4/Tue Jul 27 13:09:26 2004//
/gui_lc.m/1.4/Thu Sep 11 12:01:19 2003//
/gui_s.m/1.4/Tue Jul 27 13:09:26 2004//
/gui_sc.m/1.3/Mon Sep 8 11:28:59 2003//
/icaplot.m/1.2/Sat Apr 5 14:23:58 2003//
/pcamat.m/1.5/Mon Dec 15 18:24:32 2003//
/remmean.m/1.2/Sat Apr 5 14:23:58 2003//
/whitenv.m/1.3/Sun Oct 12 09:04:43 2003//
D
independent_component_analysis/utils/FastICA_25/Contents.m 0 → 100644
View file @ 3fb9b89e
% FastICA for Matlab 7.x and 6.x
% Version 2.5, October 19 2005
% Copyright (c) Hugo Gvert, Jarmo Hurri, Jaakko Srel, and Aapo Hyvrinen.
%
% Type fasticag to launch the graphical user interface
%
% Please refer to your Matlab documentation on how to add FastICA to your
% Matlab search path. (One place to start is the path-command)
%
% FastICA programs:
% FASTICAG - Graphical user interface for FastICA
% FASTICA - command line version of FastICA
%
% Separate functions used by FastICA programs.
% FPICA - main algorithm for calculating ICA
% WHITENV - function for whitening data
% PCAMAT - calculates the PCA for data
% REMMEAN - function for removing mean
%
% GUI_CB - needed by fasticag
% GUI_ADV - needed by fasticag
% GUI_ADVC - needed by fasticag
% GUI_L - needed by fasticag
% GUI_LC - needed by fasticag
% GUI_S - needed by fasticag
% GUI_SC - needed by fasticag
% GUI_CG - needed by fasticag
% GUI_HELP - needed by fasticag
%
% ICAPLOT - for plotting the signals
% (also used by fastica and fasticag)
%
% Misc.
% DEMOSIG - generates some test signals
%
% Deprecated
% dispsig - plots the data vectors
% replaced by icaplot
% @(#)$Id: Contents.m,v 1.7 2005/10/19 13:05:33 jarmo Exp $
\ No newline at end of file
independent_component_analysis/utils/FastICA_25/demosig.m 0 → 100644
View file @ 3fb9b89e
function [sig,mixedsig]=demosig();
%
% function [sig,mixedsig]=demosig();
%
% Returns artificially generated test signals, sig, and mixed
% signals, mixedsig. Signals are row vectors of
% matrices. Input mixedsig to FastICA to see how it works.
% @(#)$Id: demosig.m,v 1.2 2003/04/05 14:23:57 jarmo Exp $
%create source signals (independent components)
N=500; %data size
v=[0:N-1];
sig=[];
sig(1,:)=sin(v/2); %sinusoid
sig(2,:)=((rem(v,23)-11)/9).^5; %funny curve
sig(3,:)=((rem(v,27)-13)/9); %saw-tooth
sig(4,:)=((rand(1,N)<.5)*2-1).*log(rand(1,N)); %impulsive noise
for t=1:4
sig(t,:)=sig(t,:)/std(sig(t,:));
end
%remove mean (not really necessary)
[sig mean]=remmean(sig);
%create mixtures
Aorig=rand(size(sig,1));
mixedsig=(Aorig*sig);
independent_component_analysis/utils/FastICA_25/dispsig.m 0 → 100644
View file @ 3fb9b89e
function dispsig(signalMatrix, range, titlestr);
%DISPSIG - deprecated!
%
% Please use icaplot instead.
%
% See also ICAPLOT
% @(#)$Id: dispsig.m,v 1.2 2003/04/05 14:23:57 jarmo Exp $
fprintf('\nNote: DISPSIG is now deprecated! Please use ICAPLOT.\n');
if nargin < 3, titlestr = ''; end
if nargin < 2, range = 1:size(signalMatrix, 1); end
icaplot('dispsig',signalMatrix',0,range,range,titlestr);
independent_component_analysis/utils/FastICA_25/fastica.m 0 → 100644
View file @ 3fb9b89e
function [Out1, Out2, Out3] = fastica(mixedsig, varargin)
%FASTICA - Fast Independent Component Analysis
%
% FastICA for Matlab 7.x and 6.x
% Version 2.5, October 19 2005
% Copyright (c) Hugo Gävert, Jarmo Hurri, Jaakko Särelä, and Aapo Hyvärinen.
%
% FASTICA(mixedsig) estimates the independent components from given
% multidimensional signals. Each row of matrix mixedsig is one
% observed signal. FASTICA uses Hyvarinen's fixed-point algorithm,
% see http://www.cis.hut.fi/projects/ica/fastica/. Output from the
% function depends on the number output arguments:
%
% [icasig] = FASTICA (mixedsig); the rows of icasig contain the
% estimated independent components.
%
% [icasig, A, W] = FASTICA (mixedsig); outputs the estimated separating
% matrix W and the corresponding mixing matrix A.
%
% [A, W] = FASTICA (mixedsig); gives only the estimated mixing matrix
% A and the separating matrix W.
%
% Some optional arguments induce other output formats, see below.
%
% A graphical user interface for FASTICA can be launched by the
% command FASTICAG
%
% FASTICA can be called with numerous optional arguments. Optional
% arguments are given in parameter pairs, so that first argument is
% the name of the parameter and the next argument is the value for
% that parameter. Optional parameter pairs can be given in any order.
%
% OPTIONAL PARAMETERS:
%
% Parameter name Values and description
%
%======================================================================
% --Basic parameters in fixed-point algorithm:
%
% 'approach' (string) The decorrelation approach used. Can be
% symmetric ('symm'), i.e. estimate all the
% independent component in parallel, or
% deflation ('defl'), i.e. estimate independent
% component one-by-one like in projection pursuit.
% Default is 'defl'.
%
% 'numOfIC' (integer) Number of independent components to
% be estimated. Default equals the dimension of data.
%
%======================================================================
% --Choosing the nonlinearity:
%
% 'g' (string) Chooses the nonlinearity g used in
% the fixed-point algorithm. Possible values:
%
% Value of 'g': Nonlinearity used:
% 'pow3' (default) g(u)=u^3
% 'tanh' g(u)=tanh(a1*u)
% 'gauss g(u)=u*exp(-a2*u^2/2)
% 'skew' g(u)=u^2
%
% 'finetune' (string) Chooses the nonlinearity g used when
% fine-tuning. In addition to same values
% as for 'g', the possible value 'finetune' is:
% 'off' fine-tuning is disabled.
%
% 'a1' (number) Parameter a1 used when g='tanh'.
% Default is 1.
% 'a2' (number) Parameter a2 used when g='gaus'.
% Default is 1.
%
% 'mu' (number) Step size. Default is 1.
% If the value of mu is other than 1, then the
% program will use the stabilized version of the
% algorithm (see also parameter 'stabilization').
%
%
% 'stabilization' (string) Values 'on' or 'off'. Default 'off'.
% This parameter controls wether the program uses
% the stabilized version of the algorithm or
% not. If the stabilization is on, then the value
% of mu can momentarily be halved if the program
% senses that the algorithm is stuck between two
% points (this is called a stroke). Also if there
% is no convergence before half of the maximum
% number of iterations has been reached then mu
% will be halved for the rest of the rounds.
%
%======================================================================
% --Controlling convergence:
%
% 'epsilon' (number) Stopping criterion. Default is 0.0001.
%
% 'maxNumIterations' (integer) Maximum number of iterations.
% Default is 1000.
%
% 'maxFinetune' (integer) Maximum number of iterations in
% fine-tuning. Default 100.
%
% 'sampleSize' (number) [0 - 1] Percentage of samples used in
% one iteration. Samples are chosen in random.
% Default is 1 (all samples).
%
% 'initGuess' (matrix) Initial guess for A. Default is random.
% You can now do a "one more" like this:
% [ica, A, W] = fastica(mix, 'numOfIC',3);
% [ica2, A2, W2] = fastica(mix, 'initGuess', A, 'numOfIC', 4);
%
%======================================================================
% --Graphics and text output:
%
% 'verbose' (string) Either 'on' or 'off'. Default is
% 'on': report progress of algorithm in text format.
%
% 'displayMode' (string) Plot running estimates of independent
% components: 'signals', 'basis', 'filters' or
% 'off'. Default is 'off'.
%
% 'displayInterval' Number of iterations between plots.
% Default is 1 (plot after every iteration).
%
%======================================================================
% --Controlling reduction of dimension and whitening:
%
% Reduction of dimension is controlled by 'firstEig' and 'lastEig', or
% alternatively by 'interactivePCA'.
%
% 'firstEig' (integer) This and 'lastEig' specify the range for
% eigenvalues that are retained, 'firstEig' is
% the index of largest eigenvalue to be
% retained. Default is 1.
%
% 'lastEig' (integer) This is the index of the last (smallest)
% eigenvalue to be retained. Default equals the
% dimension of data.
%
% 'interactivePCA' (string) Either 'on' or 'off'. When set 'on', the
% eigenvalues are shown to the user and the
% range can be specified interactively. Default
% is 'off'. Can also be set to 'gui'. Then the user
% can use the same GUI that's in FASTICAG.
%
% If you already know the eigenvalue decomposition of the covariance
% matrix, you can avoid computing it again by giving it with the
% following options:
%
% 'pcaE' (matrix) Eigenvectors
% 'pcaD' (matrix) Eigenvalues
%
% If you already know the whitened data, you can give it directly to
% the algorithm using the following options:
%
% 'whiteSig' (matrix) Whitened signal
% 'whiteMat' (matrix) Whitening matrix
% 'dewhiteMat' (matrix) dewhitening matrix
%
% If values for all the 'whiteSig', 'whiteSig' and 'dewhiteMat' are
% supplied, they will be used in computing the ICA. PCA and whitening
% are not performed. Though 'mixedsig' is not used in the main
% algorithm it still must be entered - some values are still
% calculated from it.
%
% Performing preprocessing only is possible by the option:
%
% 'only' (string) Compute only PCA i.e. reduction of
% dimension ('pca') or only PCA plus whitening
% ('white'). Default is 'all': do ICA estimation
% as well. This option changes the output
% format accordingly. For example:
%
% [whitesig, WM, DWM] = FASTICA(mixedsig,
% 'only', 'white')
% returns the whitened signals, the whitening matrix
% (WM) and the dewhitening matrix (DWM). (See also
% WHITENV.) In FastICA the whitening matrix performs
% whitening and the reduction of dimension. Dewhitening
% matrix is the pseudoinverse of whitening matrix.
%
% [E, D] = FASTICA(mixedsig, 'only', 'pca')
% returns the eigenvector (E) and diagonal
% eigenvalue (D) matrices containing the
% selected subspaces.
%
%======================================================================
% EXAMPLES
%
% [icasig] = FASTICA (mixedsig, 'approach', 'symm', 'g', 'tanh');
% Do ICA with tanh nonlinearity and in parallel (like
% maximum likelihood estimation for supergaussian data).
%
% [icasig] = FASTICA (mixedsig, 'lastEig', 10, 'numOfIC', 3);
% Reduce dimension to 10, and estimate only 3
% independent components.
%
% [icasig] = FASTICA (mixedsig, 'verbose', 'off', 'displayMode', 'off');
% Don't output convergence reports and don't plot
% independent components.
%
%
% A graphical user interface for FASTICA can be launched by the
% command FASTICAG
%
% See also FASTICAG
% @(#)$Id: fastica.m,v 1.14 2005/10/19 13:05:34 jarmo Exp $
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Check some basic requirements of the data
if nargin == 0,
error ('You must supply the mixed data as input argument.');
end
if length (size (mixedsig)) > 2,
error ('Input data can not have more than two dimensions.');
end
if any (any (isnan (mixedsig))),
error ('Input data contains NaN''s.');
end
if ~isa (mixedsig, 'double')
fprintf ('Warning: converting input data into regular (double) precision.\n');
mixedsig = double (mixedsig);
end
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Remove the mean and check the data
[mixedsig, mixedmean] = remmean(mixedsig);
[Dim, NumOfSampl] = size(mixedsig);
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Default values for optional parameters
% All
verbose = 'on';
% Default values for 'pcamat' parameters
firstEig = 1;
lastEig = Dim;
interactivePCA = 'off';
% Default values for 'fpica' parameters
approach = 'defl';
numOfIC = Dim;
g = 'pow3';
finetune = 'off';
a1 = 1;
a2 = 1;
myy = 1;
stabilization = 'off';
epsilon = 0.0001;
maxNumIterations = 1000;
maxFinetune = 5;
initState = 'rand';
guess = 0;
sampleSize = 1;
displayMode = 'off';
displayInterval = 1;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Parameters for fastICA - i.e. this file
b_verbose = 1;
jumpPCA = 0;
jumpWhitening = 0;
only = 3;
userNumOfIC = 0;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Read the optional parameters
if (rem(length(varargin),2)==1)
error('Optional parameters should always go by pairs');
else
for i=1:2:(length(varargin)-1)
if ~ischar (varargin{i}),
error (['Unknown type of optional parameter name (parameter' ...
' names must be strings).']);
end
% change the value of parameter
switch lower (varargin{i})
case 'stabilization'
stabilization = lower (varargin{i+1});
case 'maxfinetune'
maxFinetune = varargin{i+1};
case 'samplesize'
sampleSize = varargin{i+1};
case 'verbose'
verbose = lower (varargin{i+1});
% silence this program also
if strcmp (verbose, 'off'), b_verbose = 0; end
case 'firsteig'
firstEig = varargin{i+1};
case 'lasteig'
lastEig = varargin{i+1};
case 'interactivepca'
interactivePCA = lower (varargin{i+1});
case 'approach'
approach = lower (varargin{i+1});
case 'numofic'
numOfIC = varargin{i+1};
% User has supplied new value for numOfIC.
% We'll use this information later on...
userNumOfIC = 1;
case 'g'
g = lower (varargin{i+1});
case 'finetune'
finetune = lower (varargin{i+1});
case 'a1'
a1 = varargin{i+1};
case 'a2'
a2 = varargin{i+1};
case {'mu', 'myy'}
myy = varargin{i+1};