mfbox_fasticapearson.m

toolbox for spm 5 for data, model free analysis

function [y,A,W]=mfbox_fasticapearson(mixedsig,epsilon,maxNumIterations,borderBase,borderSlope,numOfIC,firstEig,lastEig,s_verbose)
% Pearson-ICA algorithm for Independent Component Analysis
%
% This is the core program for Pearson-ICA. 
% Don't use this program directly, but 
% use the program pearson_ica that call this code.
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ 
%
% Input: 
%   mixedsig         = samples*dim matrix of the observed mixture 
%   epsilon          = scalar (convergence constant)
%   numOfIC          = integer (components to extract)
%   firstEig         = integer (largest eigenvalue to include) 
%   lastEig          = integer (smallest eigenvalue to include) 
%   maxNumIterations = integer of maximum number of iterations
%   borderBase, 
%   borderSlope      = constants for the line for switching between
%                      Pearson and tanh are:
%                      Pearson is used iff
%                      kurtosis>borderBase(1)+borderSlope(1)*skewness^2 or
%                      kurtosis<borderBase(2)+borderSlope(2)*skewness^2. 
% Output:
%   y                = samples*numOfIC matrix of the separated sources (rows)
%   A                = estimated dim*numOfIC mixing matrix
%   W                = estimated numOfIC*dim separation matrix
%
% The function estimates the independent components from given 
% multidimensional signals. Each row of the matrix mixedsig is 
% an observed signal. For optimization Hyvarinen's fixed-point algorithm,
% see http://www.cis.hut.fi/projects/ica/fastica/, is used. 
% Some of the code in this function is based on the FastICA
% package distributed under the same lisence.           
%
% example: [y, A, W]=fasticapearson(mixedsig,0.0001,100,[2.6 4],[0 1])
%
% Copyright (c) Helsinki University of Technology,
% Signal Processing Laboratory,
% Juha Karvanen, Jan Eriksson,  and Visa Koivunen.
%
% For details, see the files readme.txt
% and gpl.txt provided within the same package.
  
for i=1:5
    [y,A,W,noconv] = pearson(mixedsig,epsilon,maxNumIterations,borderBase, ...
        borderSlope,numOfIC,firstEig,lastEig,s_verbose);
    if (noconv==0), return; end
end

return


function [y,A,W,noconv]=pearson(mixedsig,epsilon,maxNumIterations,borderBase,borderSlope,numOfIC,firstEig,lastEig,s_verbose)

myaxis = gca;
Tanhconst = 1;

mixedmean = mean(mixedsig,1);
for i=1:size(mixedsig,2), mixedsig(:,i) = mixedsig(:,i)-mixedmean(i); end

[E,D,N] = mfbox_pcamat(mixedsig',firstEig,lastEig);
[whiteningMatrix,dewhiteningMatrix] = mfbox_whitenv(E,D,N);

X = real(mixedsig*whiteningMatrix');

[numSamples,Dim] = size(X);
numOfIC = min(numOfIC,Dim);

FTypeOld = zeros(numOfIC,1);
B = orth(rand(Dim,numOfIC)-.5);
BOld = B;

switch lower(s_verbose)
    case 'on'
        b_verbose = 1;
    case 'off'
        b_verbose = 0;
    otherwise
        error('Illegal value [ %s ] for parameter: verbose', s_verbose);
end

change = zeros(1,maxNumIterations+1);
score1 = zeros(size(X,1),size(B,2));
score2 = zeros(size(X,1),size(B,2));

noconv = 0;
for round=1:(maxNumIterations+1)
    if (round==(maxNumIterations+1))
        if (b_verbose>0), warning('No convergence after %d steps', maxNumIterations); end
        noconv = 1;
        break;
    end
  
    U = X*B;
    
    alpha3 = mean(U.^3);
    alpha4 = mean(U.^4);

    FType = zeros(size(U,2),1);
    for i=1:size(U,2),
        if (alpha4(i)<borderBase(1)+borderSlope(1)*alpha3(i)^2)|| ...
	        (alpha4(i)>borderBase(2)+borderSlope(2)*alpha3(i)^2),
            FType(i) = 3; %tanh
            score1(:,i) = tanh(Tanhconst*U(:,i));
            score2(:,i) = 1-score1(:,i).^2;
        else
            FType(i) = 1; %Pearson       
            [pearsona,pearsonb] = pearson_momentfit(alpha3(i),alpha4(i));
            [score1(:,i),score2(:,i)] = pearson_score(U(:,i),pearsonb,pearsona,20);
        end   
    end

    B = (X'*score1-repmat(sum(score2),size(B,1),1).*B)/numSamples;
  
    if (~all(isfinite(B(:))))
        B = randn(size(B));
        B = B*real((B'*B)^(-1/2));
    end
    B = B*real((B'*B)^(-1/2));
    while (~all(isfinite(B(:))))
        B = randn(size(B));
        B = B*real((B'*B)^(-1/2));
    end

    minAbsCos = min(abs(diag(B'*BOld)));
    change(round) = 1-minAbsCos;
    if (b_verbose>0)
        myaxis = myplot(myaxis,1:round,change(1:round));
        set(myaxis,'YScale','log'); title(myaxis,'Change per iteration');
        drawnow;
    end
           
    if ((1-minAbsCos<epsilon)&&(max(abs(FType-FTypeOld))==0)), break; end  

    BOld = B;
    FTypeOld = FType; 
end

W = B'*whiteningMatrix;
A = dewhiteningMatrix*B;
y = mixedsig*W';
m = W*mixedmean';
for i=1:size(y,2), y(:,i) = y(:,i)+m(i); end


function [phi,phi2]=pearson_score(x,b,a,C)

% PEARSON_SCORE - Calculates the score function and its derivative
% for a distribution from the Pearson system.
%
% Input:
%   x     = the realization matrix (sample values)
%   b     = the distribution parameters vector [b0 b1 b2]
%   a     = the distribution parameters vector [a0 a1]    
%   C     = the saturation limit (scalar); 
%           After saturation -C < phi(x) < C for all x
% Output:
%   phi   = the score function matrix corresponding to x  
%   phi2  = the derivative of the score function matrix corresponding to x  
%
% Copyright (c) Helsinki University of Technology,
% Signal Processing Laboratory,
% Juha Karvanen, Jan Eriksson, and Visa Koivunen.
%
% For details see the files readme.txt
% and gpl.txt provided within the same package.

% The classification is given in
%   Karvanen J. and Koivunen V.:
%   "Blind separation methods based on Pearson system and its extensions",
%   Signal Processing, 2002, to appear
%
% The estimates are given by the equations (4) and (5) in
%   Karvanen, J.,Eriksson, J., and Koivunen, V.:
%   "Pearson System Based Method for Blind Separation", 
%   Proceedings of Second International Workshop on
%   Independent Component Analysis and Blind Signal Separation,
%   Helsinki 2000, pp. 585--590

if ((b(3)~=0)&&(b(3)~=-1)),
    delta = b(2)^2/(4*b(1)*b(3));
    if (a(2)==0), ddelta = 1;
    else ddelta = 1+(delta-1)/(1+(4*b(3)/a(2)*(1+b(3)/a(2)))*delta);
    end

    if (ddelta>=1) %unbound domain
        C_1 = (-a(2)-C*b(2)-sqrt((-C*b(2)-a(2))^2+4*C*(-C*b(1)+b(2))*b(3)))/(2*C*b(3));
        C_2 = (-a(2)-C*b(2)+sqrt((-C*b(2)-a(2))^2+4*C*(-C*b(1)+b(2))*b(3)))/(2*C*b(3));
        C_3 = (a(2)-C*b(2)-sqrt((C*b(2)-a(2))^2-4*C*(C*b(1)+b(2))*b(3)))/(2*C*b(3));
        C_4 = (a(2)-C*b(2)+sqrt((C*b(2)-a(2))^2-4*C*(C*b(1)+b(2))*b(3)))/(2*C*b(3));

        if (b(2)<0) %left bound domain
            minscorelimit = max([C_1 C_2 C_3 C_4]);
            if (isreal(minscorelimit)), x = max(x,minscorelimit); end
        else %right bound domain
            maxscorelimit = min([C_1 C_2 C_3 C_4]);
            if (isreal(maxscorelimit)), x = min(x,maxscorelimit); end
        end

    elseif ((ddelta<0)||((ddelta==0)&&(b(3)>0)))
        C_1 = (-a(2)-C*b(2)-sqrt((-C*b(2)-a(2))^2+4*C*(-C*b(1)+b(2))*b(3)))/(2*C*b(3));
        C_2 = (-a(2)-C*b(2)+sqrt((-C*b(2)-a(2))^2+4*C*(-C*b(1)+b(2))*b(3)))/(2*C*b(3));
        C_3 = (a(2)-C*b(2)-sqrt((C*b(2)-a(2))^2-4*C*(C*b(1)+b(2))*b(3)))/(2*C*b(3));
        C_4 = (a(2)-C*b(2)+sqrt((C*b(2)-a(2))^2-4*C*(C*b(1)+b(2))*b(3)))/(2*C*b(3));
        CC = sort([C_1,C_2,C_3,C_4]);

        minscorelimit = CC(2);
        maxscorelimit = CC(3);

        x = min(max(x,minscorelimit),maxscorelimit);

    end
elseif ((b(3)==0)&&(b(2)~=0))  %special case: Gamma distribution
    C_1 = (b(1)*C+a(1))/(a(2)-C*b(2));
    C_2 = (-b(1)*C+a(1))/(a(2)+C*b(2));

    if (b(2)<0) %left bound Gamma distribution
        minscorelimit = max([C_1,C_2]);
        if (isreal(minscorelimit))
            x = max(x,minscorelimit); 
        end
    else %right bound Gamma distribution
        maxscorelimit = min([C_1,C_2]);
        if (isreal(maxscorelimit))
            x = min(x,maxscorelimit); 
        end
    end
end

denominator = b(1)+b(2)*x+b(3)*x.^2;
phi = -(a(2)*x-a(1))./denominator;
phi2 = -(b(1)*a(2)+b(2)*a(1)+2*a(1)*b(3)*x-a(2)*b(3)*x.^2)./(denominator.^2);


function [pearsona,pearsonb]=pearson_momentfit(alpha3,alpha4)
% PEARSON_MOMENTFIT - Estimates the parameters of the zero mean and unit
% variance distribution from the Pearson system using the third and
% forth central moments (the method of moments).
%
% Input
%   alpha3    = sample 3rd moment
%   alpha4    = sample 4th moment
%  
% Output
%   pearsona = the distribution parameters vector [a0 a1]
%   pearsonb = the distribution parameters vector [b0 b1 b2]
%
% Copyright (c) Helsinki University of Technology,
% Signal Processing Laboratory,
% Juha Karvanen, Jan Eriksson, and Visa Koivunen.
%
% For details see the files readme.txt
% and gpl.txt provided within the same package.

% The estimates are given by the equations (16)--(19) in
%   Karvanen J. and Koivunen V.:
%   "Blind separation methods based on Pearson system and its extensions",
%   Signal Processing, 2002, to appear
%
% see also
%   Karvanen, J.,Eriksson, J., and Koivunen, V.:
%   "Pearson System Based Method for Blind Separation", 
%   Proceedings of Second International Workshop on
%   Independent Component Analysis and Blind Signal Separation,
%   Helsinki 2000, pp. 585--590

if (alpha4<alpha3^2+1)
    warning('impossible values for moments %f %f',alpha3,alpha4);
end;

denominator = abs(10*alpha4-12*alpha3^2-18);
b0 = -(4*alpha4-3*alpha3^2);
b1 = -alpha3*(alpha4+3);
b2 = -(2*alpha4-3*alpha3^2-6);
pearsonb = [b0,b1,b2];
pearsona = [b1,denominator];


function a=myplot(a,varargin)

if (ishandle(a))
    if(a==gca)
        plot(a,varargin{:});
        return
    end
end
figure;
plot(varargin{:});
a = gca;