mfbox_cumindependence.m

toolbox for spm 5 for data, model free analysis

function [ci,M,A]=mfbox_cumindependence(varargin)
% measures the mutual independence of X using fourth-order cumulants
%
% usage: [ci,M,A]=mfbox_cumindependence(X[,Y],[assumewhite[,num[,normalize]]])
%
% if assumewhite is set to 1 the algo assumes white sources 
% if not (default) use C - R*E*R - trace(E*R)*R - R*E'*R; with R=Cov(X')
%
% normalization is done by dimension and by norm
%
% Copyright by Peter Gruber and Fabian J. Theis
% Signal Processing & Information Theory group
% Institute of Biophysics, University of Regensburg, Germany
% Homepage: http://research.fabian.theis.name
%           http://www-aglang.uni-regensburg.de
%
% This file is free software, subject to the 
% GNU GENERAL PUBLIC LICENSE, see gpl.txt

assumewhite = false;
normalize = false;
num = Inf;

if (nargin<1), return; end
X = varargin{1};
Y = [];
if (nargin>1)
    c = 1;
    if (length(varargin{2})==1)
        assumewhite = logical(varargin{2}); c = 2;
    else
        Y = varargin{2};
        if (nargin>2)
            if (length(varargin{3})==1)
                assumewhite = logical(varargin{3}); c = 3;
            end
        end
    end
    if (nargin>c)
        if (length(varargin{c+1})==1), num = varargin{c+1}; c = c+1; end
    end
    if (nargin>c)
        if (length(varargin{c+1})==1), normalize = logical(varargin{c+1}); end
    end
end

[T,d] = size(X); [V,e] = size(Y);

isrx = isreal(X);
if (~assumewhite),
    Xmean = mean(X,1);
    for i=1:d, X(:,i) = X(:,i)-Xmean(i); end
end

if (isempty(Y))
    if (~assumewhite)
        Rr = X'*X/(T-1);
        if (isrx), Rc = Rr;
        else, Rc = X.'*X/(T-1);
        end
    end
    isr = isrx;
else
    if (T~=V), error('incompatible data vectors'); end
    isry = isreal(Y);
    if (~assumewhite)
        Ymean = mean(Y,1);
        for i=1:e, Y(:,i) = Y(:,i)-Ymean(i); end
        Rr = [X'*X,X'*Y;Y'*X,Y'*Y]/(T-1);
        if (isry), Rc = Rr;
        else, Rc = [X.'*X,Y.'*X;X.'*Y,Y.'*Y]/(T-1);
        end
    else 
        Rr = (X'*Y)/(T-1);
        if (isry), Rc = Rr;
        else, Rc = (X'*Y)/(T-1);
        end
    end
    isr = isrx&&isry;
end

if (isempty(Y))
    M = zeros(d);
    ci = 0;
    if (nargout>2)
        if (isr), A = zeros(d,d,d*(d+1)/2);
        else A = zeros(d,d,d*d);
        end
    end
    n = 1;
    j = 0;
    for i=1:d
        if (i==j), SX = SY;
        else SX = repmat(X(:,i),1,d).*X;
        end
        if (isr), sj = [i,d:-1:(i+1)]; 
        else sj = [1,d:-1:2];
        end
        for j=sj
            % (i,j)-th crosscumulant matrix
            % should be zero in case of block independence
            if (i==j), SY = SX;
            else, SY = repmat(X(:,j),1,d).*X;
            end
            C = (SX'*SY)/T;
            if (~assumewhite)
                C = C-conj(Rc(:,i))*Rc(j,:)-Rr(:,j)*Rr(i,:)-Rr(i,j)*Rr;
                dR = sqrt(abs(diag(Rr)));
                C = C./(dR*dR'*dR(i)*dR(j));
            else
                if (i==j), C(i,j) = C(i,j)-2; C = C-eye(d);
                else, C(i,j) = C(i,j)-1; C(j,i) = C(j,i)-1;
                end
            end
            if (normalize)
                if (isr), dC = sqrt(abs(diag(C)));
                else, dC = sqrt(diag(C));
                end
                C = C./(dC*dC');
            end
            if (nargout>2), A(:,:,n) = C; end
            C = C.*conj(C);
            M = M+C; ci = ci+sum(C(logical(1-eye(d)))); % use squared frobenius norm as distance
            n = n+1;
        end
    end
    % dimension normalization
    ci = ci/n; M = M/n;
else
    M = zeros(d,e); ci = 0;
    if (nargout>2) A = zeros(d,d,d*e); end
    ui = (d+1):(d+e); li = 1:d;
    dX = ones(d,1); dY = ones(e,1);
    for i=1:d
        SX = repmat(X(:,i),1,d).*X;
        for j=1:e
            dj = d+j;
            % (i,j)-th crosscumulant matrix
            % should be zero in case of block independence
            SY = repmat(Y(:,j),1,e).*Y;
            C = (SX'*SY)/T;
            if (~assumewhite)
                C = C-conj(Rc(li,i))*Rc(dj,ui)-Rr(li,dj)*Rr(i,ui)-Rr(i,dj)*Rr(li,ui);
                dX = sqrt(abs(diag(Rr(li,li))));
                dY = sqrt(abs(diag(Rr(ui,ui))));
                C = C./(dX*dY'*dX(i)*dY(j));
            else
                C = C-conj(Rc(li,i))*Rc(dj,ui)-Rr(li,dj)*Rr(i,ui)-Rr(i,j)*Rr;
            end
            if (normalize)
                if (isr)
                    dCX = sqrt(abs(sum(X.^4,1)/T-dX'.^2));
                    dCY = sqrt(abs(sum(Y.^4,1)/T-dY'.^2));
                else
                    dCX = sqrt(sum(conj(X.^2).*(X.^2),1)/T-dX'.^2);
                    dCY = sqrt(sum(conj(Y.^2).*(Y.^2),1)/T-dY'.^2);
                end
                C = C./(dCX'*dCY);
            end
            if (nargout>2), A(:,:,((j-1)*d)+i) = C; end
            C = C.*conj(C);
            M = M+C; ci = ci+sum(C(:)); % use squared frobenius norm as distance
        end
    end
    % dimension normalization
    ci = ci/(d*e); M = M/(d*e);
end

if (nargout>2)
    if (~isr), A = conj(A); end

    s = size(A);
    if (length(s) >2 && min(num,s(3))<s(3))
        [U,D] = svd(reshape(A,[],s(3)),'econ');
        D = diag(D);
        [v,p] = sort(-abs(D));
        p = p(1:num);
        A = reshape(U(:,p)*diag(D(p)),[s(1:2),num]);
    end
end

return