# [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

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-0001.html>
```