This is a sample of python implementation for time-varying effective connectivity analysis, based on Autoregressive (AR) model with variational Bayesian noise adaptive Kalman Filter (vbKF).
\AR_kalman_connectivity
├── \figures … contains the figures generated by the script (main.py)
│
├─ \my_modules
│ ├─ ar_kalman_connectivity.py … contains main module of time-variant AR model with vbKF scheme
│ └─ myfunc.py … contains the user-defined modules
│
├── \save_data … contains the dataset of estimation results
│
└─ main.py
… contains sample script for how to use of the module “ar_kalman_connectivity.py”
Operation has been confirmed only under the following environment.
- OS: Windows 10 64bit, Ubuntu 18.04.5 LTS
- Python 3.8.3, 3.8.5
- IPython 7.16.1, 7.19.0
- conda 4.8.4, 4.9.2
- Spyder 4.1.4, 4.1.5
- numpy 1.18.5, 1.19.2
- scipy 1.5.0, 1.5.2
- matplotlib 3.2.2, 3.3.2
- scikit-learn 0.23.1, 0.23.2
- joblib 0.16.0, 0.17.0
The implemented scripts are not guaranteed to run on any other version in Python than the above.
Hiroshi Yokoyama
(Division of Neural Dynamics, Department of System Neuroscience, National Institute for Physiological Sciences, Japan)
[Simple example for use]
# Define the module
from my_modules.ar_kalman_connectivity import AR_Kalman
# Initialize
model = AR_Kalman(x, p, q, flimits)
### x .... input time-series, with size of [samples x channels]
### p .... model order for AR model
### q .... scaling factor for process noise covariance (e.g., q = 1E-5)
### flimits .... frequency band to estimate the time-variant connectivity (e.g., flimits = np.array([8, 12]) )
# Fit the model and Estmate the time-variant connectivity
model.est_kalman()
# Get the estimation results
y_hat = model.y_hat # mean of posterior predictive distribution for the observation model, with size of [samples x channels]
S = model.S # covariance of posterior predictive distribution for the observation model, with size of [channels x channels x samples]
ELBO = model.ELBO # evidence lower bound, with size of [samples x 1]
PDC = model.PDC # time-variant partial directed coherence (PDC), with size of [channels x channels x samples][Example with model selection]
# Define the module
from my_modules.ar_kalman_connectivity import AR_Kalman
import joblib
# self-defined function for use the "AR_Kalman" module with parapllel loop
def calc_model(x, p, q, flimits):
model = AR_Kalman(x, p, q, flimits)
model.est_kalman()
return model, p
# Define candidates of the parameters
P_candi = np.arange(1, 11, 1)
Qscale = np.array([10**-i for i in range(0,5)])
criteria = np.zeros((len(Qscale), len(P_candi)))
# Determine the optimal parameter
for i, q in zip(np.arange(len(Qscale)), Qscale):
#%% Calculate time-variant AR coefficients for each model order P with noise scaling factor uc
processed = joblib.Parallel(n_jobs=int(0.5*njob_max), verbose=5)(joblib.delayed(calc_model)(x, p, q, flimits) for p in P_candi)
processed.sort(key=lambda x: x[1]) # sort the output list according to the model order
tmp_result = [tmp[0] for tmp in processed]
#%% Determine the optimal model order
for j in range(len(P_candi)):
tmp = tmp_result[j]
k = tmp.Kb0.shape[0]
n = tmp.y_hat.shape[0]
elbo = tmp.ELBO.mean()# evidence lower bound
criteria[i,j] = elbo
# Select the optimal parameters
c_best = criteria.reshape(-1).max()
for i, q in zip(np.arange(len(Qscale)), Qscale):
for j, p in zip(np.arange(len(P_candi)), P_candi):
if criteria[i,j]==c_best:
Q_opt = q
P_opt = p
break
# Fit the model and Estmate the time-variant connectivity with optimal parameters
model, _ = calc_model(x, P_opt, Q_opt, flimits)
# Get the estimation results
y_hat = model.y_hat # mean of posterior predictive distribution for the observation model, with size of [samples x channels]
S = model.S # covariance of posterior predictive distribution for the observation model, with size of [channels x channels x samples]
ELBO = model.ELBO # evidence lower bound, with size of [samples x 1]
PDC = model.PDC # time-variant partial directed coherence (PDC), with size of [channels x channels x samples]Please cite my GitHub repository if you use this code in your own work:
@software{AR_kalman_connectivity,
author = {Yokoyana H.},
title = {AR_kalman_connectivity},
url = {https://github.com/myGit-YokoyamaHiroshi/AR_kalman_connectivity},
year = {2021}
}