Skip to content

Instantly share code, notes, and snippets.

@robince
Created September 20, 2023 07:51
Show Gist options
  • Save robince/de88832a931808f38d137bc1903c81c4 to your computer and use it in GitHub Desktop.
Save robince/de88832a931808f38d137bc1903c81c4 to your computer and use it in GitHub Desktop.
Bias correction permutation
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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment