omar

omar.py

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]}')