using NeuroAnalyzer
using NeuroStats
using Plots
= load("files/eeg.hdf");
eeg = epoch(eeg, ep_len=10);
e10 delete_epoch!(e10, ep=1:10);
[ Info: Loaded: EEG (24 × 308480 × 1; 1204.996 s)
Load data:
using NeuroAnalyzer
using NeuroStats
using Plots
= load("files/eeg.hdf");
eeg = epoch(eeg, ep_len=10);
e10 delete_epoch!(e10, ep=1:10);
[ Info: Loaded: EEG (24 × 308480 × 1; 1204.996 s)
(!) Most statistical functions used by the NeuroAnalyzer are in the separate package: NeuroStats.jl.
Get the first 2 seconds of the channel 1 of epochs 1:10 and calculate its 95%CI:
= e10.data[1, 1:512, 1:10]
s = e10.epoch_time[1:512]
t = bootstrap_ci(s, cl=0.95) s_avg, s_l, s_u
(s_avg = [11.654898139351815, 10.262248215819774, 11.97869012655289, 12.630062992707838, 11.570191509089227, 11.540959091716795, 10.541115219980027, 10.531517956966704, 10.421700762819713, 10.740942023467158 … -14.93467142418011, -17.482067855648076, -16.37788284377348, -12.132211640586203, -9.671487207031488, -12.30071929298087, -15.216810566070887, -15.194893417448966, -12.842532095529629, -11.009415842579445],
s_ci_l = [10.822860742097781, 9.408763518405909, 11.092509420219539, 11.649015715837571, 10.657671056665913, 10.635817998559563, 9.682447968681974, 9.764797992741789, 9.582581624381268, 9.847133243994605 … -17.115832768623285, -20.197938065265543, -19.033876224269644, -14.177775578676833, -11.463232853913574, -14.505658567550112, -17.938809501531136, -17.867950175936876, -14.920114352300565, -12.864605422159261],
s_ci_h = [12.512148422280918, 11.109775921957953, 12.871125723524765, 13.63360775095151, 12.480159298981198, 12.455200730221277, 11.419837150051526, 11.32791623654601, 11.311718997478932, 11.67202688837758 … -12.866814069530454, -14.88802574341585, -13.869626726889976, -10.211028221997651, -7.9874481181035675, -10.231731430775788, -12.669568623887448, -12.68246561037163, -10.854067971247654, -9.236188498148802],)
Plot CI:
plot_ci(s_avg, s_l, s_u, t)
Any statistical function may be calculated using bootstrapping method. The function must be provided as f="function_name(obj)"
, obj will be replaced with the signal data, e.g. f="mean(obj)"
.
= e10.data[1, :, 1:10]
s = mean(s)
s_stat = bootstrap_stat(s, f="mean(obj)") s_dist
3000-element Vector{Float64}:
0.20066477399049693
0.14394971467457118
0.07324944491156433
0.03000076203671398
0.1515606187996063
0.0641443042164969
-0.11203781864159464
0.13498272826955598
0.10262449640463948
0.05463395196986447
0.03215512833206482
0.08024724075223144
0.13264935639988237
⋮
0.026014188132643133
0.17937168246312005
0.07191929213262852
0.2452989615234543
0.2514382077903489
0.001506255250420807
-0.13471432149644222
0.2210818992910987
0.11510522773132888
0.021775359510878545
-0.04079052021458942
-0.0016443799952444494
plot_histogram(s_dist, s_stat, draw_median=false, draw_mean=true)
[ Info: Proportion of values > 0.09340081426013924: 0.529
[ Info: Proportion of values < 0.09340081426013924: 0.471