% CONDSTAT - accumulate surrogate data for comparing two data conditions
%
% Usage:
% >> [diffres, accres, res1, res2] = condstat(formula, naccu, alpha, ...
% bootside, condboot, arg1, arg2 ...);
%
% Inputs:
% formula - [string or cell array of strings] formula(s) to compute a given measure.
% Takes arguments 'arg1', 'arg2' ... and X as inputs. e.g.,
% 'sum(arg1(:,:,X),3) ./ sqrt(sum(arg2(:,:,X))) ./ sqrt(sum(arg3(:,:,X)))'
% naccu - [integer] number of accumulations of surrogate data. e.g., 200
% alpha - [float] significance level (0<alpha<0.5)
% bootside - ['both'|'upper'] side of the surrogate distribution to
% consider for significance. This parameter affect the size
% of the last dimension of accumulation array 'accres' (size
% is 2 for 'both' and 1 for 'upper').
% condboot - ['abs'|'angle'|'complex'|''] When comparing two conditions, compare
% either absolute vales ('abs'), angles ('angles') or complex values
% ('complex'). Either '' or 'complex' leave the formula unchanged;
% 'abs' takes its norm before subtraction, and 'angle' normalizes
% each value (to norm 1) before taking the difference.
% arg1 - [cell_array] of two 1D,2D,or 3D matrices of values to compare.
% The last dimension of the array is shuffled to accumulate data, the
% other dimensions must be the same size across matrices.
% e.g. size(arg1{1})=[100 200 500], size(arg1{2})=[100 200 395]
% arg2 - same as arg1, note that it is compared only to itself, and has
% nothing to do with arg1 besides using the same formula, alpha, etc.
% ...argn - may call n number of argument pairs
%
% Outputs:
% diffres - difference array for the actual (non-shuffled) data, if more than one
% arg pair is called, format is a cell array of matrices.
% accres - [cell array of 3D numerical arrays] for shuffled data, one per formula.
% res1 - result for first condition
% res2 - result for second condition
%
% Authors: Arnaud Delorme & Scott Makeig
% CNL/Salk Institute 1998-2001; SCCN/INC/UCSD, La Jolla, 2002-
%
% See also: TIMEF, CROSSF
% Copyright (C) 2002 Arnaud Delorme, Lars Kai Hansen & Scott Makeig, SCCN/INC/UCSD
%
% This file is part of EEGLAB, see http://www.eeglab.org
% for the documentation and details.
%
% Redistribution and use in source and binary forms, with or without
% modification, are permitted provided that the following conditions are met:
%
% 1. Redistributions of source code must retain the above copyright notice,
% this list of conditions and the following disclaimer.
%
% 2. Redistributions in binary form must reproduce the above copyright notice,
% this list of conditions and the following disclaimer in the documentation
% and/or other materials provided with the distribution.
%
% THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS"
% AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
% IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE
% ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT HOLDER OR CONTRIBUTORS BE
% LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
% CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
% SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
% INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
% CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
% ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF
% THE POSSIBILITY OF SUCH DAMAGE.
function [diffres, accres, res1, res2] = condstat(formula, naccu, alpha, bootside, condboot, varargin);
if nargin < 6
help condstat;
return;
end
if ~ischar(formula) && ~iscell(formula)
error('The first argument must be a string formula or cell array of string');
end
if ischar(formula)
formula = { formula };
end
if ischar(bootside)
bootside = { bootside };
end
for index = 1:length(bootside)
if ~strcmpi(bootside{index}, 'both') && ~strcmpi(bootside{index}, 'upper')
error('Bootside must be either ''both'' or ''upper''');
end
end
if ischar(condboot)
condboot = { condboot };
end
for index = 1:length(condboot)
if isempty(condboot{index}), condboot{index} = 'complex'; end
end
for index = 1:length(varargin)
if ~iscell(varargin) || length(varargin{index}) ~=2
error('Except for the first arguments, all other arguments given to the function must be cell arrays of two numerical array');
end
end
% accumulate coherence images (all arrays [nb_points x timesout x trials])
% ---------------------------
for index=1:length(varargin)
tmpvar1 = varargin{index}{1};
tmpvar2 = varargin{index}{2};
if index == 1 % Shouldn't this be recalculated for arg2, etc.? TF 2007.06.04
cond1trials = size(tmpvar1,ndims(tmpvar1));
cond2trials = size(tmpvar2,ndims(tmpvar2));
for tmpi = 1:length(formula)
accres{tmpi} = zeros(size(tmpvar1,1), size(tmpvar1,2), naccu);
end
end
if ndims(tmpvar1) == 2
eval( [ 'arg' int2str(index) '=zeros(size(tmpvar1,1), cond1trials+cond2trials);' ] );
eval( [ 'arg' int2str(index) '(:,1:cond1trials)=tmpvar1;' ] );
eval( [ 'arg' int2str(index) '(:,cond1trials+1:end)=tmpvar2;' ] );
else
eval( [ 'arg' int2str(index) '=zeros(size(tmpvar1,1), size(tmpvar1,2), cond1trials+cond2trials);' ] );
eval( [ 'arg' int2str(index) '(:,:,1:cond1trials)=tmpvar1;' ] );
eval( [ 'arg' int2str(index) '(:,:,cond1trials+1:end)=tmpvar2;' ] );
end
end
fprintf('Accumulating permutation statistics:');
alltrials = [1:cond1trials+cond2trials];
% processing formulas
% -------------------
formula1 = [];
formula2 = [];
for index = 1:length(formula)
% updating formula
% ----------------
switch lower(condboot{ min(length(condboot), index) })
case 'abs', formula{index} = [ 'abs(' formula{index} ')' ];
%case 'angle', formula{index} = [ 'exp(j*angle(' formula{index} '))' ];
case 'angle', formula{index} = [ 'angle(' formula{index} ')/(2*pi)' ];
case 'complex', ;
otherwise, condboot, error('condstat argument must be either ''abs'', ''angle'', ''complex'' or empty');
end
% computing difference (non-shuffled)
% -----------------------------------
X = 1:cond1trials;
eval( [ 'res1{index} = ' formula{index} ';'] );
X = cond1trials+1:cond1trials+cond2trials;
eval( [ 'res2{index} = ' formula{index} ';'] );
diffres{index} = res1{index}-res2{index};
% build formula to execute
% ------------------------
arrayname = [ 'accres{' int2str(index) '}' ];
if ndims(tmpvar1) == 2 % 2 dimensions
formula1 = [ formula1 arrayname '(:,index) = ' formula{index} ';'];
formula2 = [ formula2 arrayname '(:,index) = ' arrayname '(:,index) - ' formula{index} ';'];
else % 3 dimensions
formula1 = [ formula1 arrayname '(:,:,index) = ' formula{index} ';'];
formula2 = [ formula2 arrayname '(:,:,index) = ' arrayname '(:,:,index) - ' formula{index} ';'];
end
end
% accumulating (shuffling)
% -----------------------
for index=1:naccu
if rem(index,10) == 0, fprintf(' %d',index); end
if rem(index,120) == 0, fprintf('\n'); end
alltrials = shuffle(alltrials);
X = alltrials(1:cond1trials);
eval( formula1 );
X = alltrials(cond1trials+1:end);
eval( formula2 );
end
% significance level
% ------------------
for index= 1:length(formula)
accarray = accres{index};
% size = nb_points*naccu
% size = nb_points*naccu*times
if ~isreal(accarray) % might want to introduce a warning here: a complex
% result may not be desirable, and a single complex value
% in accarray could turn this from a 2-tail to 1-tail
% bootstrap test, leading to false positives. Hard to
% think of a meaningful warning though...
% TF 2007.06.04
accarray = sqrt(accarray .* conj(accarray)); % faster than abs()
end
% compute bootstrap significance level
i = round(naccu*alpha);
switch ndims(accarray)
case 3
accarray = sort(accarray,3); % always sort on naccu (when 3D, naccu is the second dim)
if strcmpi(bootside{min(length(bootside), index)}, 'upper')
accarray = accarray(:,:,floor(naccu-i/2+1));
else
accarray = accarray(:,:,[end:-1:1]);
accarraytmp(:,:,2) = accarray(:,:,ceil(i/2));
accarraytmp(:,:,1) = accarray(:,:,floor(naccu-i/2+1));
accarray = accarraytmp;
end
case 2
accarray = sort(accarray,2); % always sort on naccu (when 3D, naccu is the second dim)
if strcmpi(bootside{min(length(bootside), index)}, 'upper')
accarray = accarray(:,floor(naccu-i/2+1));
else
accarraytmp(:,2) = accarray(:,ceil(i/2));
accarraytmp(:,1) = accarray(:,floor(naccu-i/2+1));
accarray = accarraytmp;
end
case 1
accarray = sort(accarray,1); % always sort on naccu (when 3D, naccu is the second dim)
if strcmpi(bootside{min(length(bootside), index)}, 'upper')
accarray = accarray(floor(naccu-i/2+1));
else
accarraytmp(2) = accarray(ceil(i/2));
accarraytmp(1) = accarray(floor(naccu-i/2+1));
accarray = accarraytmp;
end
end
accres{index} = accarray;
end
if length(res1) == 1
res1 = res1{1};
res2 = res2{1};
diffres = diffres{1};
accres = accres{1};
end
fprintf('\n');
return;
% writing a function
% ------------------
% $$$ fid = fopen('tmpfunc.m', 'w');
% $$$ fprintf(fid, 'function [accres] = tmpfunc(alltrials, cond1trials, naccu,');
% $$$ for index=1:length(varargin)
% $$$ fprintf(fid, 'arg%d', index);
% $$$ if index ~=length(varargin), fprintf(fid,','); end
% $$$ end
% $$$ fprintf(fid, ')\n');
% $$$ commandstr = [ 'for index=1:naccu, ' ]
% $$$ % 'if rem(index,10) == 0, disp(index); end;' ];
% $$$ % 'if rem(index,10) == 0, fprintf('' %d'',index); end;' ...
% $$$ % 'if rem(index,120) == 0, fprintf(''\n''); end;' ];
% $$$ commandstr = [ commandstr 'shuffle(alltrials);' ];
% $$$ commandstr = [ commandstr 'X = alltrials(1:cond1trials);' ];
% $$$ if ndims(tmpvar1) == 2 % 2 dimensions
% $$$ commandstr = [ commandstr 'accres(:,index) = ' formula ';'];
% $$$ commandstr = [ commandstr 'X = alltrials(cond1trials+1:end);'];
% $$$ commandstr = [ commandstr 'accres(:,index) = accres(:,index)-' formula ';end;'];
% $$$ else
% $$$ commandstr = [ commandstr 'res1 = ' formula ';' 10];
% $$$ commandstr = [ commandstr 'X = alltrials(cond1trials+1:end);' 10];
% $$$ commandstr = [ commandstr 'res2 = ' formula ';' 10];
% $$$ commandstr = [ commandstr 'accres(:,:,index) = res1 - res2; end;'];
% $$$ end;
% $$$ fprintf(fid, commandstr);
% $$$ fclose(fid);
% $$$ profile on;
% $$$ accres = tmpfunc(alltrials, cond1trials, naccu, arg1);
% $$$ profile report;
% $$$ profile off;
% $$$ return;
% evaluating a command
% --------------------
% $$$ commandstr = [ 'for index=1:naccu, ' ...
% $$$ 'if rem(index,10) == 0, fprintf('' %d'',index); end;' ...
% $$$ 'if rem(index,120) == 0, fprintf(''\n''); end;' ];
% $$$ commandstr = [ commandstr 'shuffle(alltrials);' ];
% $$$ commandstr = [ commandstr 'X = alltrials(1:cond1trials);' ];
% $$$ if ndims(tmpvar1) == 2 % 2 dimensions
% $$$ commandstr = [ commandstr 'accres(:,index) = ' formula ';'];
% $$$ commandstr = [ commandstr 'X = alltrials(cond1trials+1:end);'];
% $$$ commandstr = [ commandstr 'accres(:,index) = accres(:,index)-' formula ';end;'];
% $$$ else
% $$$ commandstr = [ commandstr 'accres(:,:,index) = ' formula ';'];
% $$$ commandstr = [ commandstr 'X = alltrials(cond1trials+1:end);'];
% $$$ commandstr = [ commandstr 'accres(:,:,index) = accres(:,:,index)-' formula ';end;'];
% $$$ end;
% $$$ eval(commandstr);
% $$$ return;
% $$$
Datasets
Models
