[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