[Eeglablist] Post-AMICA utility tools (load AMICA results and various visualizations)

Jason Palmer japalmer29 at gmail.com
Tue Apr 11 06:44:51 PDT 2017


minfo.m function now pasted below (post was rejected due to attachment) …

 

From: Jason Palmer [mailto:japalmer at ucsd.edu] 
Sent: Friday, April 07, 2017 6:01 PM
To: 'Sankar Alagapan'; 'eeglablist at sccn.ucsd.edu'
Subject: RE: [Eeglablist] Post-AMICA utility tools (load AMICA results and various visualizations)

 

Hi Sankar,

 

I’m not sure if you got a response to this email, but the minfo.m file is attached.

 

Also, when you have dimensionality reduction (W size less than number of channels), the sphering matrix S is actually just the first num_pcs rows of S. So you would use, e.g.:

 

npcs = EEG.etc.amica.num_pcs;  % or npcs = size(EEG.etc.amica.W,1);

compAct = (EEG.etc.amica.W(:,:,i)*EEG.etc.amica.S(1:npcs,EEG.icachansind))*EEG.data(EEG.icachansind,:);

 

Best,

Jason

 

function [MI,vMI,MInrm,Hu,vHu,Hx,vHx] = minfo(u,x,nbins,nlags,verb)

% function MI = minfo(u,x,nbins,nlags) : Calculate the pairwise mutual

% information between rows of u, or between rows of u and rows of x if

% second argument is a matrix. If nlags > 0, calculate for lags

% -nlags:nlags. The output is a matrix of size [nu,nx,2*nlags+1] of

% pairwise mutual information values, with zero lag at MI(:,:,nlags+1).

% Uses nbins bins for pdf of u(i,:), equally spaced between min(u(i,:))

% and max(u(i,:)), similarly for pdfs of x rows.

%

% Inputs:

%           u       Matrix (nu by Nu) of nu time series.

%           x       Optional second input matrix (nx by Nx) of nx time

%                   series. Nx must equal Nu, i.e. all time series must

%                   have the same length. Enter 0 to bypass in order to

%                   set nbins or nlags.

%           nbins   Number of bins to use in computing pdfs. Default is

%                   round(3*log2(1+Nu/10)).

%           nlags   Number of lags (in addition to zero lag) to compute

%                   mutual information (before and after zero) lag.

%                   Default is 0.

%

% Outputs:

%           MI      If no second input x, or x==0, then MI is an

%                   (nu by nu by 2*nlags+1) matrix of pairwise mutual

%                   information between rows of u at lags -nlags to nlags.

%                   If both inputs u and x are present, then MI is an

%                   (nu by nx by 2*nlags+1) matrix of pairwise mutual

%                   information between rows of u and rows of x. For

%                   example, MI(i,j,nlags+1+g) is the mutual information

%                   between u(i,1+g:N+g) and u(j,1:N), or between

%                   u(i,1+g:N+g) and x(j,1:N). So a peak at g+nlags+1 in

%                   MI(i,j,:) means i lags j by g time points (i leads j by

%                   -g time points.)

%           MInrm   Matrix (nu by nu by 2*nlags+1) of "normalized" pairwise

%                   mutual information of rows of u, relative to the

%                   average of the marginal (non-differential) bin

%                   entropies. If x argument is present, normalization is

%                   always with respect to u(i,:).

%           Hu      Vector nu by 1 of differential entropies of rows of u.

%           Hx      Vector nu by 1 of differential entropies of rows of x.

%           vXX     Variance in estimate of XX

%           

varrng = 1;

if nargin < 5

    verb = 1;

end

bias_correct = 1;

secndord = 0;

[nu,Nu] = size(u);

if nargin < 3 || isempty(nbins) || nbins == 0

    nbins = round(3*log2(1+Nu/10));

end

if nargin < 4 || nlags == 0 || isempty(nlags)

    nlags = 0;

end

if nargin >= 2

    [nx,Nx] = size(x);

    if Nx == 1

       MI = zeros(nu,nu,2*nlags+1); 

        MInrm = zeros(nu,nu,2*nlags+1);

        dox = 0;

    else        

        if Nx ~= Nu

            error('input time series must be same length');

        end

        MI = zeros(nu,nx,2*nlags+1);

        MInrm = zeros(nu,nx,2*nlags+1);        

        dox = 1;

    end

else

    MI = zeros(nu,nu,2*nlags+1);

    MInrm = zeros(nu,nu,2*nlags+1);

    dox = 0;

end

 

% bin the time series

if verb

    disp('binning the time series ...'); pause(0.1);

end

Hu = zeros(nu,2*nlags+1);

deltau = zeros(nu,1);

for i = 1:nu

    if varrng

        um = mean(u(i,:));

        us = std(u(i,:));

        umax = min(max(u(i,:)),um + 5*us);

        umin = max(min(u(i,:)),um - 5*us);

    else

        umax = max(u(i,:));

        umin = min(u(i,:));

    end

    deltau(i) = (umax-umin)/nbins;

    u(i,:) = 1 + round((nbins - 1) * (u(i,:) - umin) / (umax - umin));

    u(i,:) = min(nbins,u(i,:));

    u(i,:) = max(1,u(i,:));

    

    if nlags == 0

        pmfr = diff([0 find(diff(sort(u(i,:)))) Nu])/Nu;

        Hu(i) = -sum(pmfr.*log(pmfr));

        vHu(i) = (sum(pmfr.*(log(pmfr).^2)) - Hu(i)^2) / Nu;

        if secndord

            Hu(i) = Hu(i) + (nbins-1)/(2*Nu) - (1/12/Nu^2)*(1-sum(1./pmfr));

        else

            if bias_correct

                Hu(i) = Hu(i) + (nbins-1)/(2*Nu);

            end

        end

    else

        for g = 0:nlags

            pmfr = diff([0 find(diff(sort(u(i,1+g:Nu)))) (Nu-g)])/(Nu-g);    

            Hu(i,nlags+1+g) = -sum(pmfr.*log(pmfr));

            vHu(i,nlags+1+g) = (sum(pmfr.*(log(pmfr).^2)) - Hu(i,nlags+1+g)^2) / Nu;

            if secndord

                Hu(i,nlags+1+g) = Hu(i,nlags+1+g) + (nbins-1)/(2*Nu) - (1/12/Nu^2)*(1-sum(1./pmfr));

            else

                if bias_correct

                    Hu(i,nlags+1+g) = Hu(i,nlags+1+g) + (nbins-1)/(2*Nu);

                end

            end

            if g > 0

                pmfr = diff([0 find(diff(sort(u(i,1:Nu-g)))) (Nu-g)])/(Nu-g);    

                Hu(i,nlags+1-g) = -sum(pmfr.*log(pmfr));

                vHu(i,nlags+1-g) = (sum(pmfr.*(log(pmfr).^2)) - Hu(i,nlags+1-g)^2) / Nu;

                if secndord

                    Hu(i,nlags+1-g) = Hu(i,nlags+1-g) + (nbins-1)/(2*Nu) - (1/12/Nu^2)*(1-sum(1./pmfr));

                else

                    if bias_correct

                        Hu(i,nlags+1-g) = Hu(i,nlags+1-g) + (nbins-1)/(2*Nu);

                    end

                end

            end

        end

    end

end

if dox == 1

    Hx = zeros(nx,2*nlags+1);

    deltax = zeros(nx,1);

    for i = 1:nx

        xmax = max(x(i,:));

        xmin = min(x(i,:));

        deltax(i) = (xmax-xmin)/nbins;

        x(i,:) = 1 + round((nbins - 1) * (x(i,:) - xmin) / (xmax - xmin));

 

        if nlags == 0

            pmfr = diff([0 find(diff(sort(x(i,:)))) Nx])/Nx;                        

            Hx(i) = -sum(pmfr.*log(pmfr));

            vHx(i) = (sum(pmfr.*(log(pmfr).^2)) - Hx(i)^2) / Nx;

            if secndord

                Hx(i) = Hx(i) + (nbins-1)/(2*Nx) - (1/12/Nx^2)*(1-sum(1./pmfr));

            else

                if bias_correct

                    Hx(i) = Hx(i) + (nbins-1)/(2*Nx);

                end

            end

        else

            for g = 0:nlags

                pmfr = diff([0 find(diff(sort(x(i,1+g:Nx)))) (Nx-g)])/(Nx-g);

                Hx(i,nlags+1+g) = -sum(pmfr.*log(pmfr));

                vHx(i,nlags+1+g) = (sum(pmfr.*(log(pmfr).^2)) - Hx(i,nlags+1+g)^2) / Nx;

                if secndord

                    Hx(i,nlags+1+g) = Hx(i,nlags+1+g) + (nbins-1)/(2*Nx) - (1/12/Nx^2)*(1-sum(1./pmfr));

                else

                    if bias_correct

                        Hx(i,nlags+1+g) = Hx(i,nlags+1+g) + (nbins-1)/(2*Nx);

                    end

                end

                if g > 0

                    pmfr = diff([0 find(diff(sort(x(i,1:Nx-g)))) (Nx-g)])/(Nx-g);    

                    Hx(i,nlags+1-g) = -sum(pmfr.*log(pmfr));

                    vHx(i,nlags+1-g) = (sum(pmfr.*(log(pmfr).^2)) - Hx(i,nlags+1-g)^2) / Nx;

                    if secndord

                        Hx(i,nlags+1-g) = Hx(i,nlags+1-g) + (nbins-1)/(2*Nx) - (1/12/Nx^2)*(1-sum(1./pmfr));

                    else

                        if bias_correct

                            Hx(i,nlags+1-g) = Hx(i,nlags+1-g) + (nbins-1)/(2*Nx);

                        end

                    end

                end

            end

        end

    end

else

    Hx = 0;

    vHx = 0;

end

 

% get pairwise histograms at each lag

if verb

    disp('getting pairwise mutual information ...'); pause(0.1);

end

if dox == 0

    for i = 1:nu

        if i == 1

            if verb

                disp(['doing ' int2str(i) ' ...']); pause(0.01);

            end

            tic;

        else

            t1 = toc; T = t1 * (nu-i+1);

            if verb

                disp(['doing ' int2str(i) ' ... time rem: ' num2str(T/60) ' m']); pause(0.01);

            end

            tic;

        end

        for j = 1:nu

            if nlags == 0 % faster if we don't need to compute lag index                

                if i == j

                    MI(i,i) = Hu(i);

                    MInrm(i,i) = 1;

                else

                    pmf2r = diff([0 find(diff(sort(u(i,:)+nbins*(u(j,:)-1)))) Nu])/Nu;            

                    H2 = -sum(pmf2r.*log(pmf2r));

                    if secndord

                        H2 = H2 + (nbins^2-1)/(2*Nu) - (1/12/Nu^2)*(1-sum(1./pmf2r));

                    else

                        if bias_correct

                            H2 = H2 + (nbins^2-1)/(2*Nu);

                        end                     

                    end

                    MI(i,j) = Hu(i) + Hu(j) - H2;

                    MI(j,i) = MI(i,j);

 

                    %vMI(i,j) = (sum(pmf2r.*(log(pmf2r).^2)) - H2^2) / Nu;

                    vMI(i,j) = vHu(i) + vHu(j);

                    vMI(j,i) = vMI(i,j);

 

                    MInrm(i,j) = MI(i,j) / Hu(i);

                    MInrm(j,i) = MI(j,i) / Hu(j);

                end                

            else

                % do zero lag

                if i == j

                    MI(i,i,nlags+1) = Hu(i,nlags+1);

                    MInrm(i,i,nlags+1) = 1;

                else                    

                    pmf2r = diff([0 find(diff(sort(u(i,1:Nu)+nbins*(u(j,1:Nu)-1)))) Nu])/Nu;

                    H2 = -sum(pmf2r.*log(pmf2r));

                    if secndord

                        H2 = H2 + (nbins^2-1)/(2*Nu) - (1/12/Nu^2)*(1-sum(1./pmf2r));

                    else

                        if bias_correct

                            H2 = H2 + (nbins^2-1)/(2*Nu);

                        end

                    end

                    MI(i,j,nlags+1) = Hu(i,nlags+1) + Hu(j,nlags+1) - H2;

                    MI(j,i,nlags+1) = MI(i,j,nlags+1);

 

                    %vMI(i,j,nlags+1) = (sum(pmf2r.*(log(pmf2r).^2)) - H2^2) / Nu;

                    vMI(i,j,nlags+1) = vHu(i,nlags+1) + vHu(j,nlags+1);

                    vMI(j,i,nlags+1) = vMI(i,j,nlags+1);

 

                    MInrm(i,j,nlags+1) = MI(i,j,nlags+1) / Hu(i,nlags+1);

                    MInrm(j,i,nlags+1) = MI(j,i,nlags+1) / Hu(j,nlags+1);

                end

 

                for g = 1:nlags

                    pmf2r = diff([0 find(diff(sort(u(i,1+g:Nu)+nbins*(u(j,1:(Nu-g))-1)))) Nu-g])/(Nu-g);

                    H2 = -sum(pmf2r.*log(pmf2r));

                    MI(i,j,nlags+1+g) = Hu(i,nlags+1+g) + Hu(j,nlags+1-g) - (H2 + (nbins^2-1)/2/Nu);

 

                    %vMI(i,j,nlags+1+g) = vHu(i,nlags+1+g) + vHu(j,nlags+1-g) + (sum(pmf2r.*(log(pmf2r).^2)) - H2^2) / Nu;

                    vMI(i,j,nlags+1+g) = vHu(i,nlags+1+g) + vHu(j,nlags+1-g);

                    

                    MInrm(i,j,nlags+1+g) = MI(i,j,nlags+1+g) / Hu(i,nlags+1+g);

                    if i == j

                        pmf2r = diff([0 find(diff(sort(u(i,1:Nu-g)+nbins*(u(j,1+g:Nu)-1)))) Nu-g])/(Nu-g);

                        H2 = -sum(pmf2r.*log(pmf2r));

                        if secndord

                            H2 = H2 + (nbins^2-1)/(2*Nu) - (1/12/Nu^2)*(1-sum(1./pmf2r));

                        else

                            if bias_correct

                                H2 = H2 + (nbins^2-1)/(2*Nu);

                            end

                        end

                        MI(i,j,nlags+1-g) = Hu(i,nlags+1+g) + Hu(j,nlags+1-g) - H2;

 

                        %vMI(i,j,nlags+1-g) = (sum(pmf2r.*(log(pmf2r).^2)) - H2^2) / Nu;

                        vMI(i,j,nlags+1-g) = vHu(i,nlags+1+g) + vHu(j,nlags+1-g);

                        

                        MInrm(i,j,nlags+1-g) = MI(i,j,nlags+1-g) / Hu(i,nlags+1-g);

                    else

                        MI(j,i,nlags+1-g) = MI(i,j,nlags+1+g);

                        vMI(j,i,nlags+1-g) = vMI(i,j,nlags+1+g);

                        

                        MInrm(j,i,nlags+1-g) = MI(j,i,nlags+1-g) / Hu(j,nlags+1-g);

                    end

                end                

            end        

        end

    end    

else

    for i = 1:nu

        if i == 1

            if verb

                disp(['doing ' int2str(i) ' ...']); pause(0.01);

            end

            tic;

        else

            t1 = toc; t0 = t1/(nu-i+2); T = t0 * (nu-i+2)*(nu-i+1)/2;

            if verb

                disp(['doing ' int2str(i) ' ... time rem: ' num2str(T/60) ' m']); pause(0.01);

            end

        end

        for j = 1:nx

            if nlags == 0 % faster if we don't need to compute lag index

                pmf2r = diff([0 find(diff(sort(u(i,:)+nbins*(x(j,:)-1)))) Nu])/Nu;            

                H2 = -sum(pmf2r.*log(pmf2r));

                if secndord

                    H2 = H2 + (nbins^2-1)/(2*Nu) - (1/12/Nu^2)*(1-sum(1./pmf2r));

                else

                    if bias_correct

                        H2 = H2 + (nbins^2-1)/(2*Nu);

                    end

                end

                MI(i,j) = Hu(i) + Hx(j) - H2;

 

                %vMI(i,j) = (sum(pmf2r.*(log(pmf2r).^2)) - H2^2) / Nu;

                vMI(i,j) = vHu(i) + vHx(j);

 

                MInrm(i,j) = MI(i,j) / Hu(i);

            else            

                for g = 0:nlags

                    pmf2r = diff([0 find(diff(sort(u(i,1+g:Nu)+nbins*(x(j,1:(Nu-g))-1)))) Nu-g])/(Nu-g);     

                    H2 = -sum(pmf2r.*log(pmf2r));

                    if secndord

                        H2 = H2 + (nbins^2-1)/(2*Nu) - (1/12/Nu^2)*(1-sum(1./pmf2r));

                    else

                        if bias_correct

                            H2 = H2 + (nbins^2-1)/(2*Nu);

                        end

                    end

                    MI(i,j,nlags+1+g) = Hu(i,nlags+1+g) + Hx(j,nlags+1-g) - H2;

                    

                    %vMI(i,j,nlags+1+g) = (sum(pmf2r.*(log(pmf2r).^2)) - H2^2) / Nu;

                    vMI(i,j,nlags+1+g) = vHu(i,nlags+1+g) + vHx(j,nlags+1-g);

                    

                    MInrm(i,j,nlags+1+g) = MI(i,j,nlags+1+g) / Hu(i,nlags+1+g);

                

                    if g > 0                

                        pmf2r = diff([0 find(diff(sort(u(i,1:(Nu-g))+nbins*(x(j,1+g:Nu)-1)))) Nu-g])/(Nu-g);        

                        H2 = -sum(pmf2r.*log(pmf2r));

                        if secndord

                            H2 = H2 + (nbins^2-1)/(2*Nu) - (1/12/Nu^2)*(1-sum(1./pmf2r));

                        else

                            if bias_correct

                                H2 = H2 + (nbins^2-1)/(2*Nu);

                            end

                        end

                        MI(i,j,nlags+1-g) = Hu(i,nlags+1-g) + Hx(j,nlags+1+g) - H2;

 

                        %vMI(i,j,nlags+1-g) = (sum(pmf2r.*(log(pmf2r).^2)) - H2^2) / Nu;

                        vMI(i,j,nlags+1-g) = vHu(i,nlags+1-g) + vHx(j,nlags+1+g);

                        

                        MInrm(i,j,nlags+1-g) = MI(i,j,nlags+1-g) / Hu(i,nlags+1-g);

                    end

                end

            end        

        end        

    end

end

Hu = Hu(:,nlags+1);

for i = 1:nu

    Hu(i) = Hu(i) + log(deltau(i));    

end

if dox == 1

    Hx = Hx(:,nlags+1);

    for i = 1:nx

        Hx(i) = Hx(i) + log(deltax(i));

    end

end

 

 

 

From: eeglablist-bounces at sccn.ucsd.edu [mailto:eeglablist-bounces at sccn.ucsd.edu] On Behalf Of Sankar Alagapan
Sent: Saturday, March 04, 2017 5:25 AM
To: eeglablist at sccn.ucsd.edu
Subject: Re: [Eeglablist] Post-AMICA utility tools (load AMICA results and various visualizations)

 

Hi Makoto,

 

The plugin doesn't seem to contain the code to compute mutual information (minfo). Where should we get it from?

 

Also, I run into trouble when the number of components are not the same as the number of channels.

 

compAct = (EEG.etc.amica.W(:,:,i)*EEG.etc.amica.S(:,EEG.icachansind))*EEG.data(EEG.icachansind,:);

 

My W is 114x114, S is 128x128. I had run the pop_changeweights before which took care of this discrepancy. So I just temporarily changed it to EEG.icaact. Is there a better way to do this?

 

Thank you,

Sankar

On Tue, Dec 27, 2016 at 3:00 PM <eeglablist-request at sccn.ucsd.edu> wrote:

Send eeglablist mailing list submissions to
        eeglablist at sccn.ucsd.edu

To subscribe or unsubscribe via the World Wide Web, visit
        https://sccn.ucsd.edu/mailman/listinfo/eeglablist
or, via email, send a message with subject or body 'help' to
        eeglablist-request at sccn.ucsd.edu

You can reach the person managing the list at
        eeglablist-owner at sccn.ucsd.edu

When replying, please edit your Subject line so it is more specific
than "Re: Contents of eeglablist digest..."
Today's Topics:

   1. Post-AMICA utility tools (load AMICA results and  various
      visualizations) (Makoto Miyakoshi)



---------- Forwarded message ----------
From: Makoto Miyakoshi <mmiyakoshi at ucsd.edu>
To: EEGLAB List <eeglablist at sccn.ucsd.edu>, Jason Palmer <japalmer29 at gmail.com>
Cc: 
Bcc: 
Date: Mon, 26 Dec 2016 15:37:04 -0800
Subject: [Eeglablist] Post-AMICA utility tools (load AMICA results and various visualizations)

Dear lab and list,

 

Upon Scott's request, I collected AMICA utility functions that, I guess, Ozgur Balkan developed mostly. This contains AMICA result loader and various visualizations, including pairwise mutual information (PMI) plotter with diagonalization function. It can be downloaded from EEGLAB's plugin manager.

 

-- 

Makoto Miyakoshi
Swartz Center for Computational Neuroscience
Institute for Neural Computation, University of California San Diego

_______________________________________________
eeglablist mailing list eeglablist at sccn.ucsd.edu
Eeglablist page: http://www.sccn.ucsd.edu/eeglab/eeglabmail.html
To unsubscribe, send an empty email to eeglablist-unsub at sccn.ucsd.edu
To switch to non-digest mode, send an empty email to eeglablist-nodigest at sccn.ucsd.edu

-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://sccn.ucsd.edu/pipermail/eeglablist/attachments/20170411/dcfcb162/attachment.html>


More information about the eeglablist mailing list