subj = 15; % CHANGE THIS NUMBER FOR EACH SUBJECT

%% load the raw data from disk
prefix = sprintf('Sub%02d', subj);
load([prefix '_raw']);
load([prefix '_avg_Faces_vs_Scrambled']);

%% reorganize the timelocked data and compute stats

cfg = [];
cfg.channel    = 'MEGMAG';
cfg.keeptrials = 'yes';
timelock = ft_timelockanalysis(cfg, data);

cfg           = [];
cfg.correctm  = 'no';
cfg.method    = 'analytic';
cfg.statistic = 'indepsamplesT';        % this is implemented in ft_statfun_indepsamplesT
cfg.design    = nan(1, size(timelock.trialinfo,1));
cfg.design(timelock.trialinfo==1) = 1;  % Famous faces
cfg.design(timelock.trialinfo==2) = 1;  % Unfamiliar faces
cfg.design(timelock.trialinfo==3) = 2;  % Scrambled
cfg.ivar      = 1;                      % the first (and only) row of the design represents the independent variable
analytic = ft_timelockstatistics(cfg, timelock);


%% do some sanity checks
figure
imagesc(analytic.time, 1:length(analytic.label), -log10(analytic.prob))
colorbar

cfg = [];
cfg.channel = analytic.label;
tmp = ft_selectdata(cfg, avg_Faces_vs_Scrambled);

analytic.avg = tmp.avg;

% analytic.logprob = -log10(analytic.prob);
% analytic.logprob(isnan(analytic.logprob)) = 0;
% analytic.logprob(isinf(analytic.logprob)) = 10;

% save analytic analytic

cfg               = [];
cfg.layout        = 'neuromag306mag.lay';
cfg.parameter     = 'avg';
cfg.maskparameter = 'mask';
figure; ft_multiplotER(cfg, analytic);

%% use montecarlo and correctm=max

cfg                  = [];
cfg.correctm         = 'max';
cfg.method           = 'montecarlo';
cfg.numrandomization = 1000;
cfg.statistic        = 'indepsamplesT'; % this is implemented in ft_statfun_indepsamplesT
cfg.design           = nan(1, size(timelock.trialinfo,1));
cfg.design(timelock.trialinfo==1) = 1; % Famous faces
cfg.design(timelock.trialinfo==2) = 1; % Unfamiliar faces
cfg.design(timelock.trialinfo==3) = 2; % Scrambled
cfg.ivar             = 1; % the first (and only) row of the design represents the independent variable
cfg.latency          = [0.140 0.180];
montecarlo = ft_timelockstatistics(cfg, timelock);

% save montecarlo montecarlo

figure
hist([montecarlo.negdistribution' montecarlo.posdistribution'], 100)

%% compare the observed statistical values to the distributions

negdistribution = sort(montecarlo.negdistribution);
negthreshold    = negdistribution(26)   % why not at 5%, i.e. 51

posdistribution = sort(montecarlo.posdistribution);
posthreshold    = posdistribution(975)  % why not at 5%, i.e. 950

figure
subplot(3,1,1)
hist(montecarlo.negdistribution, 50)
ylabel('negdist');
xlim([-10 10]);
subplot(3,1,2)
hist(montecarlo.stat(:), 100)
ylabel('observed stat');
xlim([-10 10]);
subplot(3,1,3)
hist(montecarlo.posdistribution, 50)
ylabel('posdist');
xlim([-10 10]);


%% use your own trialfunction, e.g. spearman rank correlation

cfg                  = [];
cfg.channel          = 'MEG2021';
cfg.statistic        = 'statfun_parametric';
cfg.design           = nan(1, size(timelock.trialinfo,1));
cfg.design(timelock.trialinfo==1) = 1; % Famous faces
cfg.design(timelock.trialinfo==2) = 2; % Unfamiliar faces
cfg.design(timelock.trialinfo==3) = 3; % Scrambled
cfg.ivar             = 1; % the first (and only) row of the design represents the independent variable

cfg.correctm         = 'no'; % or another method
cfg.method           = 'analytic';
analytic2 = ft_timelockstatistics(cfg, timelock);

cfg.correctm         = 'max';
cfg.method           = 'montecarlo';
cfg.numrandomization = 1000;
montecarlo2 = ft_timelockstatistics(cfg, timelock);

figure
hold on
plot(analytic2.time,   -log10(analytic2.prob),   'b')
plot(montecarlo2.time, -log10(montecarlo2.prob), 'r')
line([montecarlo2.time(1) montecarlo2.time(end)], [1.3 1.3])

% save analytic2 analytic2
% save montecarlo2 montecarlo2
