estimate_dynamics_model

 Estimate linear dynaics model constrained by anatomical connectivity
 from ROI-current estimated from MEG data

 Copyright (C) 2011, ATR All Rights Reserved.
 License : New BSD License(see VBMEG_LICENSE.txt)

0001 function estimate_dynamics_model(p)
0002 % Estimate linear dynaics model constrained by anatomical connectivity
0003 % from ROI-current estimated from MEG data
0004 %
0005 % Copyright (C) 2011, ATR All Rights Reserved.
0006 % License : New BSD License(see VBMEG_LICENSE.txt)
0007 
0008 disp(mfilename);
0009 
0010 % Set parameters
0011 lamda=1e-4;% Regularization parameter
0012 order_auto = 2;% Order of autoregressive term
0013 
0014 % Make directory to save dynamics model
0015 save_dir = fullfile(p.proj_root, p.dynamics_dirname, p.dynamics_model_dirname);
0016 if exist(save_dir, 'dir') ~= 7
0017     mkdir(save_dir);
0018 end
0019 
0020 % Load delays between ROIs
0021 dmri_file = fullfile(p.proj_root, p.dmri_dirname, p.connectivity_dirname, p.connectivity_file);
0022 load(dmri_file, 'delay_ms_matrix')
0023 
0024 for ta = 1:length(p.task_list)
0025     
0026     % Load ROI-current
0027     current_file = fullfile(p.proj_root, p.dynamics_dirname, p.roi_current_dirname, [p.task_list{ta} '.mat']);
0028     load(current_file, 'Z', 'sampling_rate')
0029     
0030     % Make delay matrix
0031     Delta=make_Delta(delay_ms_matrix, sampling_rate);
0032     
0033     % Estimate linear dynamics model constrained by anatomical connectivity
0034     while 1
0035         
0036         % Estimate linear dynamics model
0037         [A, B] = lcd_fitl2reg_holocalAR(Z, Delta, order_auto, lamda*max(sum(Z.^2,2)));
0038         
0039         % Check whether divergence occurs
0040         conMAR = lcd_ab2conmar(A, B, Delta, 0);
0041         omega = lcd_calc_residuals(Z, conMAR);
0042         if max(abs(Z(:)))<max(abs(omega(:)))
0043             lamda=lamda*10;
0044             disp(['Try Lamda=' num2str(lamda)])
0045         else
0046             disp(['Lamda is determined to be ' num2str(lamda)])
0047             break
0048         end
0049     end
0050     
0051     % Save estimated dynamics model
0052     save(fullfile(save_dir, [p.task_list{ta} '.lcd.mat']),'A', 'B', 'Delta')
0053 end
0054