Контекст и объяснение: Ниже приведен сам код, без данных, которые у меня есть для него. Я пытаюсь заставить for-l oop работать вокруг события причинности Грейнджера для некоторых данных, так как он будет запускать один и тот же анализ для пяти различных переменных. Он принимает два входа и обеспечивает один ключевой вывод, целиком код. Я пытаюсь сделать это для какой-то работы, которую я запланировал в будущем, в основном как побочный проект, но я не могу во всем мире выяснить, где разместить сингл для l oop, чтобы при построении моих переменных я мог подготовьте их все вместе и сравните их оттуда. Благодарим за любую идею! В данный момент я просто перебираю все циклы, которые мне нужны для генерации сравнений, которые, я думаю, могут быть сделаны.
function [TimeFrequency_GC] = GC_20200405(multi_channel, srate)
%IMPORTANT: At this point in time, this function was only built to take in
%EEG time series data.
% -----Function Summary------
% Granger Causality is the statistical hypothesis test that assesses whether past
%and present values of a set X time series variables referred to as the
%"cause" variables, affect the predictive distribution of a distinct set of
%Y time series variables referred to as the "effect" variables.
% For users new to statistics, this means that your null hypothesis is: X
%does not granger cause Y.
% Within the code, a dataset, variables, and loops and commands are used to
%perform the granger causality analysis.
%Mock EEG data was used to run this function.
%References
%[1] Ding, M., Chen, Y., & Bressler, S. L. (2006). 17 Granger causality:
% basic theory and application to neuroscience. Handbook of time series
% analysis: recent theoretical developments and applications, 437.
%
%[2] Seth, A. K., Barrett, A. B., & Barnett, L. (2015). Granger causality
% analysis in neuroscience and neuroimaging. Journal of Neuroscience,
% 35(8), 3293-3297.
%
%[3] Timme, N. M., & Lapish, C. (2018).
% A tutorial for information theory in neuroscience. ENeuro, 5(3).
%%
%Initialize variables and paramters
% Set channels
A = 'O1';
B = 'PZ';
C = 'P3';
D = 'P4';
E = 'FZ';
%%
%Use a laplacian transformation to attain time derivatives *Taken from the
%book
EEG.data = laplacian_perrinX(EEG.data, [EEG.chanlocs.X],[EEG.chanlocs.Y],[EEG.chanlocs.Z]);
%The time window for the data bins and the model order is measured in
%milliseconds, as well as the downsampling rate. Downsampling (decimation)
%helps in reducing data size, compression, or image reduction. This helps
%to conserve memory and signal processing time when the data have been
%oversampled (which may often be the case), thus improving the results.
TimeWindow = 100; ModelOrder = 10; Decimation = -400:20:1200;
srate = 1000/EEG.srate; %Here is the sampling rate, in milliseconds
%Next, the indices are converted into parameters
TimeWindowPoints = round(TimeWindow/(srate));
Model_Order_Points = round(ModelOrder/(srate));
%Locate the index variables within the EEG data set
A = find(strcmpi(A,{EEG.chanlocs.labels}));
B = find(strcmpi(B,{EEG.chanlocs.labels}));
C = find(strcmpi(C,{EEG.chanlocs.labels}));
D = find(strcmpi(D,{EEG.chanlocs.labels}));
E = find(strcmpi(E,{EEG.chanlocs.labels}));
multi_channel = [A B C D E];
%Without maintaining stationarity, the analysis loses accuracy. Therefore,
%event related potentials are removed that disrupt stationarity.
avgEEGdata = mean(EEG.data((multi_channel),:,:), 3);
EEGdata = bsxfun(@minus,EEG.data(multi_channel,:,:), avgEEGdata);
Time_Indices = dsearchn(EEG.times', Decimation');
%Here, the variables selected are compared to each other, in both ways.
%Then, the Bayes Info Criterion is calculated.
[a2b,b2a] = deal(zeros(1, length(Decimation)));
BayesInfoCriterion = zeros(length(Decimation), ModelOrder);
%%
%For-loop #1: This loop will create a time index for each decimation
%(downsampled data point)
% the for-loop I want should go here, I assume, and I thought of setting it up as for i = 1:length(multi_channel). Next, I thought that anywhere I have another indexing variable I should then put (i) right after it to loop through all of the channels.
for Timeindex = 1:length(Decimation)
TimeFrequency_GC = squeeze(EEGdata(:, Time_Indices(Timeindex)- floor(TimeWindowPoints/2): Time_Indices(Timeindex) + floor(TimeWindowPoints/2) - mod(TimeWindowPoints+1,2),:));
%For-loop #2: This loop will Z-score and detrend all data points
for TrialIndex = 1:size(TimeFrequency_GC)
TimeFrequency_GC(1, :, TrialIndex) = zscore(detrend(squeeze(TimeFrequency_GC(1, :, TrialIndex))));
end
%Next, the data are reshaped to fit the armorf function borrowed from
%the book's code.
EEGTimePoints = TimeWindowPoints*EEG.trials;
TimeFrequency_GC = reshape(TimeFrequency_GC, 5, EEGTimePoints);
% This line was taken from the book for simplicity
[Aa,Ex] = armorf(TimeFrequency_GC(1,:), EEG.trials, TimeWindowPoints, Model_Order_Points);
[Ab,Ey] = armorf(TimeFrequency_GC(2,:), EEG.trials, TimeWindowPoints, Model_Order_Points);
[Ac,Ex] = armorf(TimeFrequency_GC(3,:), EEG.trials, TimeWindowPoints, Model_Order_Points);
[Ad,Ey] = armorf(TimeFrequency_GC(4,:), EEG.trials, TimeWindowPoints, Model_Order_Points);
[Ae,Ey] = armorf(TimeFrequency_GC(5,:), EEG.trials, TimeWindowPoints, Model_Order_Points);
[Aab,E] = armorf(TimeFrequency_GC, EEG.trials, TimeWindowPoints, Model_Order_Points);
%Next, calculate the time-domain causal estimate
%The following lines were modified to represent multiple channels from
%one to the next.
a2b(Timeindex)=log(Ey/E(2,2));
a2c(Timeindex)=log(Ey/E(3,3));
a2d(Timeindex)=log(Ey/E(4,4));
a2e(Timeindex)=log(Ey/E(5,5));
%The following lines represent the reveresed flow of information within
%specific channels.
b2a(Timeindex)=log(Ex/E(1,1));
c2a(Timeindex)=log(Ey/E(2,2));
d2a(Timeindex)=log(Ey/E(3,3));
e2a(Timeindex)=log(Ey/E(4,4));
%
% A Bayes Information Criteria (BIC) is generated here to locate optimal
% model orders at each time point within the dataset. A BIC is often
% used for neural datasets as it compensates for the large number of
% data points found within neuroelectric datasets.
%For-loop #3: This loop will identify Bayes Information Criteria for
%optimal model functionality. In turn, the user can (A) select the optimal
%model order in their data set or (B) average the model orders to attempt
%to generate a well-rounded model order. NOTE: With Bayesian/Schwarz
%Information Criterions, the lowest criterion is often preferred, as it
%reflects the lowest difference between all of the models generated.
for BayesInfoCriteriaIndex = 1:size(BayesInfoCriterion,5)
% Run the model and compute the Bayes Information Criterion
[Aab,E] = armorf(TimeFrequency_GC,EEG.trials,TimeWindowPoints,BayesInfoCriteriaIndex);
BayesInfoCriterion(Timeindex,BayesInfoCriteriaIndex) = log(det(E)) ...
+ (log(length(TimeFrequency_GC))*BayesInfoCriteriaIndex*2^2)/length(TimeFrequency_GC);
%
end
end
%Initialize plot variables for robust plotting
BayesInfo = (1:size(BayesInfoCriterion,5));
[BayesInfoCriteriaVal,BayesInfoCriteriaIdx] = min(mean(BayesInfoCriterion, 1));
BestOrder = (BayesInfoCriteriaIdx*(1000/EEG.srate));
%Plot figures, with labels, titles, and legend(s) for comparison and
%analysis
figure(1)
plot(Decimation, a2b)
hold on
plot(Decimation, a2c)
plot(Decimation, a2d)
plot(Decimation, a2e)
title('Multiple Forward Granger Causality Analysis')
xlabel('Time (ms)')
ylabel('Granger Estimation Values')
figure(2)
plot(Decimation, b2a)
hold on
plot(Decimation, c2a)
plot(Decimation, d2a)
plot(Decimation, e2a)
title('Multiple Reverse Granger Causality Analysis')
xlabel('Time (ms)')
ylabel('Granger Estimation Values')
%This second subplot will show the Bayes Info Criterion, from its first
%point to the size of itself.
figure(3)
plot((BayesInfo)*(srate), mean(BayesInfoCriterion, 2),'--.')
xlabel('Order (ms)')
ylabel('Average Bayes Info Criterion for all time points')
hold on
plot(BayesInfoCriteriaIdx*(srate), BayesInfoCriteriaVal, 'k*','markersize',12)
title('Best Model Order Number for Analysis')
end