Bias correction permutation

gistfile1.txt

Ntrl_class = 720;
% can change effect strength here
effect_mean = 0.2;
% effect_mean = 0.1;
r = rng(999);
% stimulus values (ie 0: standard, 1: deviant)
y = [zeros(Ntrl_class,1); ones(Ntrl_class,1)];
% class 0 data: standard norml
x0 = randn(Ntrl_class,1);
% class 1 data: a small difference in mena
x1 = randn(Ntrl_class,1)+effect_mean;
x = [x0;x1];

% calculate MI values from the data with and without bias correction
cx = copnorm(x);
Ibc = mi_model_gd(cx,y,2,true)
Inobc = mi_model_gd(cx,y,2,false)

% permutation test
% shuffle y to break relationship with x and obtain samples
% from a surrogate empirical null distribution with the same 
% properties as the data
Nperm = 1000;
Iperm_bc = zeros(Nperm,1);
Iperm_nobc = zeros(Nperm,1);
for pi=1:Nperm
    % shuffle y to break relationship
    idx = randperm(length(y));
    yshf = y(idx);
    Iperm_bc(pi) = mi_model_gd(cx,yshf,2,true);
    Iperm_nobc(pi) = mi_model_gd(cx,yshf,2,false);
end

% calculate the percentile of the MI value measured from intact data
% with respect to the permutation surrogate null distribution
[f,x] = ecdf(Iperm_bc);
[~,m] = min(abs(x-Ibc));
Ibc_prc = f(m);
[f,x] = ecdf(Iperm_nobc);
[~,m] = min(abs(x-Inobc));
Inobc_prc = f(m);

% plot results
figure
ax = [];
ax(1) = subplot(2,1,1);
hist(Iperm_bc)
xline(Ibc,'r')
title(sprintf('with BC, percentile: %d', round(100*Ibc_prc)));
ax(2) = subplot(2,1,2);
hist(Iperm_nobc)
xline(Inobc,'r')
title(sprintf('no BC, percentile: %d', round(100*Inobc_prc)));

linkaxes(ax)