Commit deeebbac authored by Francesco Rea's avatar Francesco Rea
Browse files

adding the m-files in the proper repo

parents
function ccaStartup
%---------------------------------------------------------------
% ccaStartup.m
% add matlab paths for causal connectivity analysis toolbox
% run before anything else!
% AKS Apr 09 2009
%---------------------------------------------------------------
addpath test/;
addpath utilities/;
warning('off','MATLAB:dispatcher:InexactMatch');
% check all files are included
if exist('armorf.m') ~=2 | exist('pwcausal.m') ~=2,
error('Functions armorf.m and/or pwcausal.m missing. Download from www.brain-smart.org');
end
% check that mex file is compiled properly
dummy = 1;
try
x = marsaglia_normcdf(dummy);
catch
error('Mex file marsaglia_normcdf is not properly compiled for this operating system. Please recompile by choosing a compiler (mex -setup) and then entering: mex utilities/marsaglia_normcdf.c');
end
if abs(x-0.8413)>1e-2,
abs(x-0.8413)
error('Mex file marsaglia_normcdf is not properly compiled for this operating system (runs but gives wrong value). Please recompile by choosing a compiler (mex -setup) and then entering: mex utilities/marsaglia_normcdf.c');
end
function [ret] = cca_autonomy_regress(X,nlags)
% -----------------------------------------------------------------------
% FUNCTION: cca_autonomy_regress.m
% PURPOSE: perform multivariate regression with granger causality statistics
%
% INPUT: X - nvar (rows) by nobs (cols) observation matrix
% NLAGS - number of lags to include in model
%
% OUTPUT: ret.prb - probability with which variable i is autonomous
% w.r.t. the other variables in the set
% ret.fs - F-statistics for the same
% ret.gaut - log ratio based magnitude of autonomy
%
% Modified AKS February 2007 to implement the autonomy test
% Modified AKS December 2008 to reorganize output and do ratio test
% Modified AKS August 2009 to check residual autocorrelation &
% consistency
%
% Ref: Seth, A.K. (in press). Measuring autonomy and emergence via
% Granger causality. Artificial Lie
% Ref: Seth, A.K., Measuring autonomy via multivariate autoregressive
% modelling, in Proceedings of the 9th European Conference on Artificial Life,
% 2007, p. 475-484.
% COPYRIGHT NOTICE AT BOTTOM
% -----------------------------------------------------------------------
% figure regression parameters
nobs = size(X,2);
nvar = size(X,1);
if(nvar>nobs) error('error in cca_autonomy_regress: nvar>nobs, check input matrix'); end
% remove sample means if present (no constant terms in this regression)
X = cca_rm_temporalmean(X);
disp('computing G-autonomies');
% construct lag matrices
lags = -999*ones(nvar,nobs-nlags,nlags);
for jj=1:nvar
for ii=1:nlags
lags(jj,:,nlags-ii+1) = X(jj,ii:nobs-nlags+ii-1);
end
end
% unrestricted regression (no constant term)
regressors = zeros(nobs-nlags,nvar*nlags);
for ii=1:nvar
s1 = (ii-1)*nlags+1;
regressors(:,s1:s1+nlags-1) = squeeze(lags(ii,:,:));
end
for ii=1:nvar
xvec = X(ii,:)';
xdep = xvec(nlags+1:end);
beta(:,ii) = regressors\xdep;
xpred(:,ii) = regressors*beta(:,ii); % keep hold of predicted values
u(:,ii) = xdep-xpred(:,ii);
RSS1(ii) = sum(u(:,ii).^2);
C(ii) = covariance(u(:,ii),u(:,ii),nobs-nlags);
end
covu = cov(u);
% A rectangular matrix A is rank deficient if it does not have linearly independent columns.
% If A is rank deficient, the least squares solution to AX = B is not unique.
% The backslash operator, A\B, issues a warning if A is rank deficient and
% produces a least squares solution that has at most rank(A) nonzeros.
% restricted regressions (no constant terms)
for ii=1:nvar
xvec = X(ii,:)';
xdep = xvec(nlags+1:end); % dependent variable
eq_inx = setdiff(1:nvar,ii); % vars to include in restricted regression
regressors = zeros(nobs-nlags,length(eq_inx)*nlags);
for kk=1:length(eq_inx)
s1 = (kk-1)*nlags+1;
regressors(:,s1:s1+nlags-1) = squeeze(lags(eq_inx(kk),:,:));
end
beta_r = regressors\xdep;
u_r(:,ii) = xdep-regressors*beta_r;
S(ii) = covariance(u_r(:,ii),u_r(:,ii),nobs-nlags); % dec 08
RSS0(ii) = sum(u_r(:,ii).^2);
end
covr = cov(u_r);
% do Granger f-tests
prb = ones(nvar,1).*NaN;
ftest = zeros(nvar,1).*NaN;
gaut = ones(nvar,1).*NaN;
n2 = (nobs-nlags)-(nvar*nlags);
for ii=1:nvar
ftest(ii) = ((RSS0(ii)-RSS1(ii))/nlags)/(RSS1(ii)/n2);
prb(ii) = 1 - cca_cdff(ftest(ii),nlags,n2);
gaut(ii) = log(S(ii)/C(ii));
end
% do r-squared and check whiteness
df_error = (nobs-nlags)-(nvar*nlags);
df_total = (nobs-nlags);
for ii = 1:nvar
xvec = X(ii,nlags+1:end);
rss2 = xvec*xvec';
rss(ii) = 1 - (RSS1(ii) ./ rss2);
rss_adj(ii) = 1 - ((RSS1(ii)/df_error) / (rss2/df_total) );
waut(ii) = cca_whiteness(X,u(:,ii));
end
% check consistency
cons = cca_consistency(X,xpred,nvar);
% organize output structure
ret.gaut = single(gaut');
ret.fs = single(ftest');
ret.prb = single(prb');
ret.covu = covu;
ret.covr = covr;
ret.rss = rss;
ret.rss_adj = rss_adj;
ret.waut = waut;
ret.cons = cons;
% This file is part of GCCAtoolbox.
% It is Copyright (C) Anil Seth (2007-2009)
%
% GCCAtoolbox is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% GCCAtoolbox is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GCCAtoolbox. If not, see <http://www.gnu.org/licenses/>.
% find best organization of subplot panels
function [nx,ny] = calc_panels(nvar)
if nvar == 2,
nx = 1;
ny = 2;
elseif nvar < 5,
nx = 2;
ny = 2;
elseif nvar < 10,
nx = 3;
ny = 3;
elseif nvar < 16,
nx = 4;
ny = 4;
else
t = sqrt(nvar);
nx = ceil(t);
ny = ceil(t);
end
\ No newline at end of file
function [ret] = cca_causaldensity(GC,PR)
%-----------------------------------------------------------------------
% FUNCTION: cca_causaldensity.m
% PURPOSE: calculate causal density given a matrix of Granger causalities
%
% INPUTS: GC: nvar by nvar matrix of Granger magnitudes
% PR: nvar by nvar binary matrix of significant links (optional)
% (if PR is present, GC will be filtered by PR)
%
% OUTPUT: ret.cd: causal density (bounded in range 0:1)
% ret.cdw: causal density weighted by magnitude (not bounded)
% ret.ucd: vector of causal densities (per node)
% ret.ucdw: as above but weighted
%
% Written by Anil Seth, March 2004
% Updated August 2004 AKS
%
% Ref: Seth, A.K. (2005) Network: Comp. Neural. Sys. 16(1):35-55
%
% updated Apr 24 2006 AKS, to ensure that all sig connections = 1
% updated May 28 2008 AKS, for compatability
% updated April 2009 AKS, add ucd/ucdw outputs, reorganize I/O
% structure
% COPYRIGHT NOTICE AT BOTTOM
%-----------------------------------------------------------------------
nvar = length(GC);
if nargin == 1, % just GC matrix
ret.cd = -1;
ret.cdw = sum(sum(GC))/(nvar*(nvar-1));
ret.ucd = -1;
ret.ucdw = (sum(GC)/nvar) + ((sum(GC'))/nvar);
else % filter by significance
PR(PR~=0) = 1; % ensure that all sig connections have value = 1
ret.cd = sum(sum(PR))/(nvar*(nvar-1));
ret.ucd = (sum(PR)/nvar) + ((sum(PR'))/nvar);
GC(PR==0)=0;
ret.cdw = sum(sum(GC))/(nvar*(nvar-1));
ret.ucdw = (sum(GC)/nvar) + ((sum(GC'))/nvar);
end
% This file is part of GCCAtoolbox.
% It is Copyright (C) Anil Seth (2005-2009)
%
% GCCAtoolbox is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% GCCAtoolbox is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GCCAtoolbox. If not, see <http://www.gnu.org/licenses/>.
function [ret] = cca_causaldensity_spectral(GW,thresh)
%-----------------------------------------------------------------------
% FUNCTION: cca_causaldensity_spectral.m
% PURPOSE: calculate causal density given a matrix of Granger causalities
%
% INPUTS: GW: nvar by nvar by nfreq matrix of Granger magnitudes
% thresh: EITHER 1*nfreq vector of thresholds for GC
% OR nvar*nvar*nfreq binary matrix showing significances
% (derived via bootstrap or permutation)
%
% OUTPUT:
% ret.scdw: causal density weighted by magnitude (not
% bounded, 1*nfreq)
% ret.sucdw: unit causal density, nvar * nfreq
%
% ret.scd: 1*nfreq vector of causal densities, unweighted
% ret.sucd: unit causal density, nvar * nfreq
% (the latter two only returned for bstrap/permute matrices)
%
% AKS April 2009 AKS, based on cca_causaldensity.m
% COPYRIGHT NOTICE AT BOTTOM
%-----------------------------------------------------------------------
% check inputs
nvar = size(GW,1);
nfreq = size(GW,3);
if nargin == 1,
error('insufficient inputs to cca_causaldensity_spectral.m');
end
if min(size(thresh)) == 1,
METHOD = 1;
disp('cca_causaldensity_spectral: applying vector of thresholds, 1*nfreq');
else
METHOD = 2;
disp('cca_causaldensity_spectral: applying bootstrap/permutation significance matrix');
end
if METHOD == 1,
scdw = zeros(1,nfreq);
sucdw = zeros(nvar,nfreq);
for i=1:nfreq,
X = GW(:,:,i);
X(X<thresh(i)) = 0;
scdw(i) = sum(sum(X))/(nvar*(nvar-1));
sucdw(:,i) = ((sum(X)/nvar) + ((sum(X'))/nvar))';
end
ret.scdw = scdw;
ret.sucdw = sucdw;
else
scdw = zeros(1,nfreq);
sucdw = zeros(nvar,nfreq);
scd = zeros(1,nfreq);
sucd = zeros(nvar,nfreq);
PR = thresh;
clear thresh;
PR(PR~=0) = 1; % ensure that all sig connections have value = 1
for i=1:nfreq,
X = GW(:,:,i);
P = PR(:,:,i);
scd(i) = sum(sum(P))/(nvar*(nvar-1));
sucd(i) = (sum(P)/nvar) + ((sum(P'))/nvar);
X(P==0) = 0;
scdw(i) = sum(sum(X))/(nvar*(nvar-1));
sucdw(:,i) = ((sum(X)/nvar) + ((sum(X'))/nvar))';
end
ret.scdw = scdw;
ret.sucdw = sucdw;
ret.scd = scd;
ret.sucd = sucd;
end
% This file is part of GCCAtoolbox.
% It is Copyright (C) Anil Seth (2009)
%
% GCCAtoolbox is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% GCCAtoolbox is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GCCAtoolbox. If not, see <http://www.gnu.org/licenses/>.
function [ret] = cca_causalflow(GC,PR)
%-----------------------------------------------------------------------
% FUNCTION: cca_causalflow.m
% PURPOSE: calculate causal flow given a matrix of Granger causalities
%
% INPUTS: GC: nvar by nvar matrix of Granger magnitudes
% PR: nvar by nvar binary matrix of significant links (optional)
% (if PR is present, GC will be filtered by PR)
%
% OUTPUT: ret.indeg: only incoming causal influences
% ret.outdeg: only outgoing causal influences
% ret.flow: outdeg minus indeg
% ret.windeg: as above but weighted
% ret.woutdeg: as above but weighted
% ret.wflow: as above but weighted
%
% Written by Anil K Seth, March 2004
% Updated AKS August 2004
% Updated AKS December 2005
% Updated AKS May 2008 for compatability
% Updated AKS Apr 2009 for new toolbox, no longer uses F stat
%
% Ref: Seth, A.K. (2005) Network: Comp. Neural. Sys. 16(1):35-55
% COPYRIGHT NOTICE AT BOTTOM
%-----------------------------------------------------------------------
if nargin == 1, % just the GC matrix
ret.windeg = sum(GC');
ret.woutdeg = sum(GC);
ret.wflow = ret.woutdeg-ret.windeg;
else
PR(PR~=0) = 1; % ensure that all sig connections have value = 1 (Apr 24 2006)
% unweighted
ret.indeg = sum(PR');
ret.outdeg = sum(PR);
ret.flow = ret.outdeg-ret.indeg;
% weighted
GC(PR==0) = 0;
ret.windeg = sum(GC');
ret.woutdeg = sum(GC);
ret.wflow = ret.woutdeg-ret.windeg;
end
% This file is part of GCCAtoolbox.
% It is Copyright (C) Anil Seth (2004-2009)
%
% GCCAtoolbox is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% GCCAtoolbox is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GCCAtoolbox. If not, see <http://www.gnu.org/licenses/>.
function [ret] = cca_causalflow_spectral(GW,thresh)
%-----------------------------------------------------------------------
% FUNCTION: cca_causalflow_spectral.m
% PURPOSE: calculate causal flow given a matrix of Granger causalities,
% modified for spectral measures
%
% INPUTS: GW: nvar by nvar by nfreq matrix of Granger magnitudes
% thresh: EITHER 1*nfreq vector of thresholds for GC
% OR nvar*nvar*nfreq binary matrix showing significances
% (derived via bootstrap or permutation)
%
% OUTPUT:
% ret.swindeg: weighted incoming causal influences
% ret.swoutdeg: weighted outgoing causal influences
% ret.swflow: weighted out-in flow
%
% ret.sindeg: weighted incoming causal influences
% ret.woutdeg: weighted outgoing causal influences
% ret.sflow: weighted out-in flow
% (last 3 only if bstrap/permute matrix is provided)
%
% AKS April 2009,
% based on aks_causalflow.m
% COPYRIGHT NOTICE AT BOTTOM
%-----------------------------------------------------------------------
% check inputs
nvar = size(GW,1);
nfreq = size(GW,3);
if nargin == 1,
error('insufficient inputs to aks_causalflow_spectral.m');
end
if min(size(thresh)) == 1,
METHOD = 1;
disp('aks_causalflow_spectral: applying vector of thresholds, 1*nfreq');
else
METHOD = 2;
disp('aks_causalflow_spectral: applying bootstrap/permutation significance matrix');
end
if METHOD == 1,
windeg = zeros(nvar,nfreq);
woutdeg = zeros(nvar,nfreq);
wflow = zeros(nvar,nfreq);
for i=1:nfreq,
X = GW(:,:,i);
X(X<thresh(i)) = 0;
windeg(:,i) = sum(X');
woutdeg(:,i) = sum(X);
wflow(:,i) = woutdeg(:,i) - windeg(:,i);
end
ret.swindeg = windeg;
ret.swoutdeg = woutdeg;
ret.swflow = wflow;
else
windeg = zeros(nvar,nfreq);
woutdeg = zeros(nvar,nfreq);
wflow = zeros(nvar,nfreq);
indeg = zeros(nvar,nfreq);
outdeg = zeros(nvar,nfreq);
flow = zeros(nvar,nfreq);
PR = thresh;
clear thresh;
PR(PR~=0) = 1; % ensure that all sig connections have value = 1
for i=1:nfreq,
X = GW(:,:,i);
P = PR(:,:,i);
indeg(:,i) = sum(P');
outdeg(:,i) = sum(P);
flow(:,i) = outdeg(:,i) - indeg(:,i);
X(P==0) = 0;
windeg(:,i) = sum(X');
woutdeg(:,i) = sum(X);
wflow(:,i) = woutdeg(:,i) - windeg(:,i);
end
ret.swindeg = windeg;
ret.swoutdeg = woutdeg;
ret.swflow = wflow;
ret.sindeg = indeg;
ret.soutdeg = outdeg;
ret.sflow = flow;
end
% This file is part of GCCAtoolbox.
% It is Copyright (C) Anil Seth (2007-2009)
%
% GCCAtoolbox is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% GCCAtoolbox is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GCCAtoolbox. If not, see <http://www.gnu.org/licenses/>.
function unit_root = cca_check_cov_stat(X,NLAGS)
%-----------------------------------------------------------------------
% FUNCTION: cca_check_cov_stat.m
% PURPOSE: verify that all time-series in a multivariate set are
% covariance-stationary. Uses the Dickey-Fuller test for
% unit-roots
%
% INPUTS: X: matrix of nvar variables by nobs observations of each variable
% NLAGS: number of lags over which to make estimations (<<nobs)
%
% OUTPUT: unit_root: vector of length nvar containing '0' for no unit
% roots, and 1 for a unit root. Presence of a unit root
% indicates the corresponding time-series is not
% covariance-stationary.
%
% Written by Anil K Seth, March 2004
% Updated August 2004
% Ref: Seth, A.K. (2005) Network: Comp. Neural. Sys. 16(1):35-55
% COPYRIGHT NOTICE AT BOTTOM
%-----------------------------------------------------------------------
[nvar,nobs] = size(X);
if(nobs<nvar) error('Fewer observations than variables, exiting'); end
if(nobs<NLAGS) error('NLAGS too large, exiting'); end
unit_root = zeros(nvar,1);
for i=1:nvar
y = X(i,:);
for nlags = 1:NLAGS
res = cca_adf(y',0,nlags); % augmented Dickey-Fuller test
tstat(i,nlags) = res.adf;
cval(i,nlags) = res.crit(3); % use 10% quintile, see cca_adf.m for more info
if(tstat(i,nlags) > cval(i,nlags)) unit_root(i) = 1; end
end
end
% This file is part of GCCAtoolbox.
% It is Copyright (C) Anil Seth (2004-2009)
%
% GCCAtoolbox is free software: you can redistribute it and/or modify
% it under the terms of the GNU General Public License as published by
% the Free Software Foundation, either version 3 of the License, or
% (at your option) any later version.
%
% GCCAtoolbox is distributed in the hope that it will be useful,
% but WITHOUT ANY WARRANTY; without even the implied warranty of
% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
% GNU General Public License for more details.
%
% You should have received a copy of the GNU General Public License
% along with GCCAtoolbox. If not, see <http://www.gnu.org/licenses/>.
function unit_root = cca_check_cov_stat_mtrial(X,Nr,Nl,NLAGS)
%-----------------------------------------------------------------------
% FUNCTION: cca_check_cov_stat_mtrial.m
% PURPOSE: verify that all time-series in a multivariate set are
% covariance-stationary. Uses the Dickey-Fuller test for
% unit-roots. Adapted for multiple trials
%
% INPUTS: X: matrix of nvar variables by nobs observations of each variable
% Nr: number of realizations
% Nl: number of data points in each realization
% NLAGS: number of lags over which to make estimations (<Nr)
%
% OUTPUT: unit_root: matrix of size nvar*Nr containing '0' for no unit