Skip to content

Instantly share code, notes, and snippets.

@kingjr
Created January 3, 2022 15:11
Show Gist options
  • Save kingjr/4a7f6491bca240ed2683c54e6524501f to your computer and use it in GitHub Desktop.
Save kingjr/4a7f6491bca240ed2683c54e6524501f to your computer and use it in GitHub Desktop.
omar
import numpy as np
from scipy.stats import ttest_1samp, wilcoxon
def ttest(a):
return ttest_1samp(a, 0)
n_times = 30
n_chans = 4
n_epochs = 200
n_subjects = 100
def make_subject():
epo = np.random.randn(n_epochs, n_chans, n_times)
wfreq = np.random.rand(n_epochs)
alpha = np.random.rand(n_epochs)
return epo, wfreq, alpha
def gfp(epo):
n_trials, n_chans, n_times = epo.shape
evo = epo.mean(0)
return evo.std(0)
out = np.empty((n_subjects, n_times, 4))
for subject in range(n_subjects):
epo, wfreq, alpha = make_subject()
# factorize
wfreq_bin = wfreq>np.median(wfreq)
alpha_bin = alpha>np.median(alpha)
out[subject, :, 0] = gfp(epo[~alpha_bin & ~wfreq_bin])
out[subject, :, 1] = gfp(epo[alpha_bin & ~wfreq_bin])
out[subject, :, 2] = gfp(epo[~alpha_bin & wfreq_bin])
out[subject, :, 3] = gfp(epo[alpha_bin & wfreq_bin])
windows = [slice(0, 10), slice(10, 20), slice(20, 30)]
for stat in (wilcoxon, ttest):
print(stat.__name__)
for window in windows:
# wfreq low
alphalow = out[:, window, 0].mean(1)
alphahigh = out[:, window, 1].mean(1)
wfreqlow = alphalow - alphahigh
print(f'wfreq low: alpha low vs alpha high: p={stat(wfreqlow)[1]}')
# wfreq high
alphalow = out[:, window, 2].mean(1)
alphahigh = out[:, window, 3].mean(1)
wfreqhigh = alphalow - alphahigh
print(f'wfreq high: alpha low vs alpha high: p={stat(wfreqhigh)[1]}')
print(f'interactaion wfreq * alpha p={stat(wfreqhigh - wfreqlow)[1]}')
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment