efica.m

ICALAB for Signal Processing Toolbox for BSS, ICA, cICA, ICA-R, SCA and MCA

    case 4 %'pow3'
        mu=mean(s.^4,2);
        nu=3*ones(dim,1);
        beta=mean(s.^6,2);
end
J=ones(1,dim); gam=(nu-mu)'; jm=(beta-mu.^2)';
Werr=(jm'*J+J'*jm+J'*gam.^2)./(abs(gam)'*J+J'*abs(gam)).^2;
Werr=Werr-diag(diag(Werr));
ISRsymm=Werr/N;
%SIRsymm=-10*log10(sum(Werr,2)'/N);

ekurt=mean(pwr(s,4),2); %%% estimate fourth moment

%Werr=Werr/N+diag((ekurt-1)/(4*N));
% GCov=diag(Werr(:));
% for i=1:dim-1
%     for j=i+1:dim
%         Gjiij=(-jm(i)-jm(j)+abs(gam(i)*gam(j)))/(abs(gam(i))+abs(gam(j)))^2;
%         GCov((i-1)*dim+j,(j-1)*dim+i)=Gjiij;
%         GCov((j-1)*dim+i,(i-1)*dim+j)=Gjiij;
%     end
% end
%varW_symm=reshape(diag(kron(Wsymm',eye(dim))*GCov*kron(Wsymm,eye(dim))),dim,dim);

%varA_symm=reshape(diag(kron(eye(dim),inv(Wsymm))*GCov*kron(eye(dim),inv(Wsymm)')),dim,dim);

%EFICA


% if eficaGaussianNL=='npnl'
%     ekurt=ones(1,dim)*4;
% end
for j=1:dim
    w=W(:,j);
    if ekurt(j)<2.4184   %%% sub-Gaussian -> try "GGD score function" alpha=>3
        if ekurt(j)<1.7  %%% the distribution seems to be extremly subgaussian
            %alpha=50; % bimodal-gaussian distribution will be considered
            alpha=15;
        else %%% estimate shape parameter alpha from the fourth moment
            if ekurt(j)>1.8
                alpha=1/(sqrt((5*ekurt(j)-9)/6/pi^2)+(1.202)*3*(5*ekurt(j)-9)/pi^4);
            else
                alpha=15; %the distribution is likely uniform
            end
        end
        if alpha<3 %%% the distribution is rather gaussian -> keep the original result
            w=W_before_decor(:,j);
        elseif alpha<=15 %%% try score function sign(x)*|x|^(alpha-1)
            wold=zeros(dim,1);
            nit=0;
            alpha=ceil(min([alpha 15]));
            while ((1-abs(w'*wold/norm(w))>fineepsilon) && (nit<EFICAMaxIt) &&...
                    (abs((W(:,j)/norm(W(:,j)))'*(w/norm(w)))>min_correlation))
                w=w/norm(w);
                %abs((W(:,j)/norm(W(:,j)))'*(w/norm(w)))
                wold=w;
                u=Z'*w;
                ualpha=pwr(abs(u),alpha-2);
                w=Z*(u.*ualpha)/N-(alpha-1)*mean(ualpha)*w;
                nit=nit+1;
            end
            sest=(w/norm(w))'*Z;
            if abs((W(:,j)/norm(W(:,j)))'*(w/norm(w)))>min_correlation,
                mu(j)=mean(pwr(abs(sest),alpha));
                nu(j)=(alpha-1)*mean(pwr(abs(sest),alpha-2));
                beta(j)=mean(pwr(abs(sest),2*alpha-2));
            end
        else  %trying bimodal distribution
            wold=zeros(dim,1);
            nit=0;
            while ((1-abs(w'*wold/norm(w))>fineepsilon) && (nit<EFICAMaxIt) &&...
                    (abs((W(:,j)/norm(W(:,j)))'*(w/norm(w)))>min_correlation))
                w=w/norm(w);
                wold=w;
                u=Z'*w;
                m=mean(abs(u)); %estimate the distance of distribution's maximas
                e=sqrt(1-m^2); %then their variance is..., because the overall variance is 1
                if e<=0.05, e=0.05; m=sqrt(1-e^2); end %due to stability
                uplus=u+m; uminus=u-m;
                expplus=exp(-uplus.^2/2/e^2);
                expminus=exp(-uminus.^2/2/e^2);
                expb=exp(-(u.^2+m^2)/e^2);
                g=-(uminus.*expminus + uplus.*expplus)./(expplus+expminus)/e^2;
                gprime=-(e^2*(expplus.^2+expminus.^2)+(2*e^2-4*m^2)*expb)./(expplus+expminus).^2/e^4;
                w=Z*g/N-mean(gprime)*w;
                nit=nit+1;
            end
            u=(w/norm(w))'*Z;
            if abs((W(:,j)/norm(W(:,j)))'*(w/norm(w)))>min_correlation,
                m=mean(abs(u)); e=sqrt(1-m^2);
                uplus=u+m; uminus=u-m;
                expplus=exp(-uplus.^2/2/e^2); expminus=exp(-uminus.^2/2/e^2); expb=exp(-(u.^2+m^2)/e^2);
                g=-(uminus.*expminus + uplus.*expplus)./(expplus+expminus)/e^2;
                gprime=-(e^2*(expplus.^2+expminus.^2)+(2*e^2-4*m^2)*expb)./(expplus+expminus).^2/e^4;
                mu(j)=mean(u.*g);
                nu(j)=mean(gprime);
                beta(j)=mean(g.^2);
            end
        end %bi-modal variant
    elseif ekurt(j)>3.762      %%% supergaussian (alpha<1.5)
        switch eficaGaussianNL
            case 1 %'npnl' %nonparametric density model; the same as in NPICA algorithm; very slow
                wold=zeros(dim,1);
                nit=0;
                while ((1-abs(w'*wold/norm(w))>fineepsilon) && (nit<EFICAMaxIt) &&...
                        (abs((W(:,j)/norm(W(:,j)))'*(w/norm(w)))>min_correlation))
                    w=w/norm(w);
                    [G GP]=nonln(w'*Z);
                    wold=w;
                    w=Z*G'/N-mean(GP)*w;
                    nit=nit+1;
                end
                [G GP]=nonln((w/norm(w))'*Z);
                mu(j)=mean(((w/norm(w))'*Z).*G);nu(j)=mean(GP);beta(j)=mean(G.^2);
            case 2 %'gaus'
                wold=zeros(dim,1);
                nit=0;
                while ((1-abs(w'*wold/norm(w))>fineepsilon) && (nit<EFICAMaxIt) &&...
                        (abs((W(:,j)/norm(W(:,j)))'*(w/norm(w)))>min_correlation))
                    w=w/norm(w);
                    aexp=exp(-(w'*Z).^2/2);
                    wold=w;
                    w=Z*((w'*Z).*aexp)'/N-w*mean((1-(w'*Z).^2).*aexp);
                    nit=nit+1;
                end
                sest=(w/norm(w))'*Z; aexp=exp(-sest.^2/2);
                mu(j)=mean(sest.^2.*aexp); nu(j)=mean((1-sest.^2).*aexp);
                beta(j)=mean((sest.*aexp).^2);
            case 3 %'ggda' %our proposal
                gam=3.3476;
                wold=zeros(dim,1);
                nit=0;
                while ((1-abs(w'*wold/norm(w))>fineepsilon) && (nit<EFICAMaxIt) &&...
                        (abs((W(:,j)/norm(W(:,j)))'*(w/norm(w)))>min_correlation))
                    w=w/norm(w);
                    u=w'*Z;
                    ualpha=abs(u);
                    aexp=exp(-gam*ualpha);
                    wold=w;
                    w=Z*(u.*aexp)'/N-w*mean((1-gam*ualpha).*aexp);
                    nit=nit+1;
                end
                if abs((W(:,j)/norm(W(:,j)))'*(w/norm(w)))>min_correlation,
                    sest=(w/norm(w))'*Z; ualpha=abs(sest); aexp=exp(-gam*ualpha);
                    mu(j)=mean(sest.^2.*aexp); nu(j)=mean((1-gam*ualpha).*aexp);
                    beta(j)=mean((sest.*aexp).^2);
                end
        end
    else  % alpha \in (1.5, 3) ; keep the original nonlinearity
        w=W_before_decor(:,j);
    end
    if abs((W(:,j)/norm(W(:,j)))'*(w/norm(w)))<min_correlation
        %%%%  the signal changed too much, thus,
        %%%%  seems to be rather gaussian -> keep the original
        %%%%  nonlinearity and result
        W(:,j)=W_before_decor(:,j);
    else
        W(:,j)=w;
    end
end

%Refinement
J=abs(mu-nu)';
B=(beta-mu.^2)';
%EJ=abs(Smu-Snu)';
%EB=(Sbeta-Smu.^2)';
Werr=zeros(dim); Hef=zeros(dim);
for k=1:dim
    ccc=(J*B(k)/J(k))./(B+J.^2);ccc(k)=1;
    Hef(k,:)=(ccc.*J-ccc(k)*J(k))./(ccc.*J+ccc(k)*J(k));
    WW=W*diag(ccc);
    sloupce=setdiff(1:dim,find(sum(WW.^2)/max(sum(WW.^2))<1e-7)); % removing almost zero rows
    M=WW(:,sloupce);
    M=M*real(inv(M'*M)^(1/2));
    if sum(sloupce==k)==1
        Wefica(:,k)=M(:,sloupce==k);
    else
        w=null(M');
        if size(w,2)==1
            Wefica(:,k)=w(:,1);
        else % there are probably more than two gaussian components => use the old result to get a regular matrix
            Wefica(:,k)=Wsymm(:,k);
        end
    end
    %Estimate variance of elements of the gain matrix
    Werr(k,:)=B(k)*(B+J.^2)./(J.^2*B(k)+J(k)^2*(B+J.^2));
end
Werr=Werr-diag(diag(Werr));
ISRef=Werr/N;
%SIRefica=-10*log10(sum(Werr,2)'/N);

% Werr=Werr/N+diag((ekurt-1)/(4*N));
% GCov=diag(Werr(:));
% for i=1:dim-1
%     for j=i+1:dim
%         Gjiij=-B(i)*B(j)/(J(i)^2*(B(j)+J(j)^2)+J(j)^2*B(i));
%         GCov((i-1)*dim+j,(j-1)*dim+i)=Gjiij;
%         GCov((j-1)*dim+i,(i-1)*dim+j)=Gjiij;
%     end
% end
Wefica=Wefica'*CC;

%varW_efica=reshape(diag(kron(Wefica',eye(dim))*GCov*kron(Wefica,eye(dim))),dim,dim);

%varA_efica=reshape(diag(kron(eye(dim),inv(Wefica))*GCov*kron(eye(dim),inv(Wefica'))),dim,dim);

function [g,gprime]=nonln(X)
[d N]=size(X);
g=zeros(d,N);
gprime=zeros(d,N);
h=1.06*N^(-1/5);
for i=1:d
    for k=1:N
        Z=(X(i,k)-X(i,:))/h;
        Zsquare=Z.^2;
        Phi=exp(-Zsquare/2)/sqrt(2*pi);
        f=mean(Phi,2)'/h;
        fprime=-mean(Z.*Phi,2)'/h^2;
        fprime2=mean((Zsquare-1).*Phi,2)'/h^3;
        g(i,k)=-fprime./f;
        gprime(i,k)=-fprime2./f+g(i,k)^2;
    end
end

function x=pwr(a,n)
x=a;
for i=2:n
    x=x.*a;
end