Homec4science
Diffusion Slepians (master) Browse

Slepi/Code/Utilities/spm12/spm_dcm_estimate.m
8a186cd82567master

Last Change
Raw File

spm_dcm_estimate.m
View Options

function [DCM] = spm_dcm_estimate(P)
% Estimates parameters of a DCM (bilinear or nonlinear) for fMRI data
% FORMAT [DCM] = spm_dcm_estimate(DCM)
% DCM - DCM structure or its filename
%
% Expects
%--------------------------------------------------------------------------
% DCM.a % switch on endogenous connections
% DCM.b % switch on bilinear modulations
% DCM.c % switch on exogenous connections
% DCM.d % switch on nonlinear modulations
% DCM.U % exogenous inputs
% DCM.Y.y % responses
% DCM.Y.X0 % confounds
% DCM.Y.Q % array of precision components
% DCM.n % number of regions
% DCM.v % number of scans
%
% Options
%--------------------------------------------------------------------------
% DCM.options.two_state % two regional populations (E and I)
% DCM.options.stochastic % fluctuations on hidden states
% DCM.options.centre % mean-centre inputs
% DCM.options.nonlinear % interactions among hidden states
% DCM.options.nograph % graphical display
% DCM.options.induced % switch for CSD data features
% DCM.options.P % starting estimates for parameters
% DCM.options.hidden % indices of hidden regions
% DCM.options.nmax % maximum number of (effective) nodes
% DCM.options.nN % maximum number of iterations
%
% Evaluates:
%--------------------------------------------------------------------------
% DCM.M % Model structure
% DCM.Ep % Condition means (parameter structure)
% DCM.Cp % Conditional covariances
% DCM.Vp % Conditional variances
% DCM.Pp % Conditional probabilities
% DCM.H1 % 1st order hemodynamic kernels
% DCM.H2 % 2nd order hemodynamic kernels
% DCM.K1 % 1st order neuronal kernels
% DCM.K2 % 2nd order neuronal kernels
% DCM.R % residuals
% DCM.y % predicted data
% DCM.T % Threshold for Posterior inference
% DCM.Ce % Error variance for each region
% DCM.F % Free-energy bound on log evidence
% DCM.ID % Data ID
% DCM.AIC % Akaike Information criterion
% DCM.BIC % Bayesian Information criterion
%
%__________________________________________________________________________
% Copyright (C) 2002-2012 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_dcm_estimate.m 6432 2015-05-09 12:58:12Z karl $
SVNid = '$Rev: 6432 $';
%-Load DCM structure
%--------------------------------------------------------------------------
if ~nargin
%-Display model details
%----------------------------------------------------------------------
Finter = spm_figure('GetWin','Interactive');
set(Finter,'name','Dynamic Causal Modelling')
%-Get DCM
%----------------------------------------------------------------------
[P, sts] = spm_select(1,'^DCM.*\.mat$','select DCM_???.mat');
if ~sts, DCM = []; return; end
spm('Pointer','Watch')
spm('FigName','Estimation in progress');
end
if isstruct(P)
DCM = P;
P = ['DCM-' date '.mat'];
else
load(P)
end
% check options
%==========================================================================
try, DCM.options.two_state; catch, DCM.options.two_state = 0; end
try, DCM.options.stochastic; catch, DCM.options.stochastic = 0; end
try, DCM.options.nonlinear; catch, DCM.options.nonlinear = 0; end
try, DCM.options.centre; catch, DCM.options.centre = 0; end
try, DCM.options.nmax; catch, DCM.options.nmax = 8; end
try, DCM.options.nN; catch, DCM.options.nN = 32; end
try, DCM.options.hidden; catch, DCM.options.hidden = []; end
try, DCM.n; catch, DCM.n = size(DCM.a,1); end
try, DCM.v; catch, DCM.v = size(DCM.Y.y,1); end
try, M.nograph = DCM.options.nograph; catch, M.nograph = spm('CmdLine');end
try, M.Nmax = DCM.options.nN ; catch,M.Nmax = 64; end
try, M.Nmax = DCM.M.Nmax ; catch,M.Nmax = 64; end
% analysis and options
%--------------------------------------------------------------------------
DCM.options.induced = 0;
% unpack
%--------------------------------------------------------------------------
U = DCM.U; % exogenous inputs
Y = DCM.Y; % responses
n = DCM.n; % number of regions
v = DCM.v; % number of scans
% detrend outputs (and inputs)
%--------------------------------------------------------------------------
Y.y = spm_detrend(Y.y);
if DCM.options.centre
U.u = spm_detrend(U.u);
end
% check scaling of Y (enforcing a maximum change of 4%
%--------------------------------------------------------------------------
scale = max(max((Y.y))) - min(min((Y.y)));
scale = 4/max(scale,4);
Y.y = Y.y*scale;
Y.scale = scale;
% check confounds (add constant if necessary)
%--------------------------------------------------------------------------
if ~isfield(Y,'X0'),Y.X0 = ones(v,1); end
if ~size(Y.X0,2), Y.X0 = ones(v,1); end
% fMRI slice time sampling
%--------------------------------------------------------------------------
try, M.delays = DCM.delays; catch, M.delays = ones(n,1); end
try, M.TE = DCM.TE; end
% create priors
%==========================================================================
% check DCM.d (for nonlinear DCMs)
%--------------------------------------------------------------------------
try
DCM.options.nonlinear = logical(size(DCM.d,3));
catch
DCM.d = zeros(n,n,0);
DCM.options.nonlinear = 0;
end
% specify parameters for spm_int_D (ensuring updates every second or so)
%--------------------------------------------------------------------------
if DCM.options.nonlinear
M.IS = 'spm_int_D';
M.nsteps = round(max(Y.dt,1));
M.states = 1:n;
else
M.IS = 'spm_int';
end
% check for endogenous DCMs, with no exogenous driving effects
%--------------------------------------------------------------------------
if isempty(DCM.c) || isempty(U.u)
DCM.c = zeros(n,1);
DCM.b = zeros(n,n,1);
U.u = zeros(v,1);
U.name = {'null'};
end
if ~any(spm_vec(U.u)) || ~any(spm_vec(DCM.c))
DCM.options.stochastic = 1;
end
% priors (and initial states)
%--------------------------------------------------------------------------
[pE,pC,x] = spm_dcm_fmri_priors(DCM.a,DCM.b,DCM.c,DCM.d,DCM.options);
str = 'Using specified priors\n';
try, M.P = DCM.options.P; end % initial parameters
try, pE = DCM.options.pE; fprintf(str); end % prior expectation
try, pC = DCM.options.pC; fprintf(str); end % prior covariance
try, M.P = DCM.M.P; end % initial parameters
try, pE = DCM.M.pE; fprintf(str); end % prior expectation
try, pC = DCM.M.pC; fprintf(str); end % prior covariance
% eigenvector constraints on pC for large models
%--------------------------------------------------------------------------
if n > DCM.options.nmax
% remove confounds and find principal (nmax) modes
%----------------------------------------------------------------------
y = Y.y - Y.X0*(pinv(Y.X0)*Y.y);
V = spm_svd(y');
V = V(:,1:DCM.options.nmax);
% remove minor modes from priors on A
%----------------------------------------------------------------------
j = 1:(n*n);
V = kron(V*V',V*V');
pC(j,j) = V*pC(j,j)*V';
end
% hyperpriors over precision - expectation and covariance
%--------------------------------------------------------------------------
hE = sparse(n,1) + 6;
hC = speye(n,n)/128;
i = DCM.options.hidden;
hE(i) = -4;
hC(i,i) = exp(-16);
% complete model specification
%--------------------------------------------------------------------------
M.f = 'spm_fx_fmri'; % equations of motion
M.g = 'spm_gx_fmri'; % observation equation
M.x = x; % initial condition (states)
M.pE = pE; % prior expectation (parameters)
M.pC = pC; % prior covariance (parameters)
M.hE = hE; % prior expectation (precisions)
M.hC = hC; % prior covariance (precisions)
M.m = size(U.u,2);
M.n = size(x(:),1);
M.l = size(x,1);
M.N = 32;
M.dt = 16/M.N;
M.ns = v;
% nonlinear system identification (nlsi)
%==========================================================================
if ~DCM.options.stochastic
% nonlinear system identification (Variational EM) - deterministic DCM
%----------------------------------------------------------------------
[Ep,Cp,Eh,F] = spm_nlsi_GN(M,U,Y);
% predicted responses (y) and residuals (R)
%----------------------------------------------------------------------
y = feval(M.IS,Ep,M,U);
R = Y.y - y;
R = R - Y.X0*spm_inv(Y.X0'*Y.X0)*(Y.X0'*R);
Ce = exp(-Eh);
else
% proceed to stochastic (initialising with deterministic estimates)
%======================================================================
% Decimate U.u from micro-time
% ---------------------------------------------------------------------
u = U.u;
y = Y.y;
Dy = spm_dctmtx(size(y,1),size(y,1));
Du = spm_dctmtx(size(u,1),size(y,1));
Dy = Dy*sqrt(size(y,1)/size(u,1));
u = Dy*(Du'*u);
% DEM Structure: place model, data, input and confounds in DEM
% ---------------------------------------------------------------------
DEM.M = M;
DEM.Y = y';
DEM.U = u';
DEM.X = Y.X0';
% set inversion parameters
% ---------------------------------------------------------------------
DEM.M(1).E.form = 'Gaussian'; % form of random fluctuations
DEM.M(1).E.s = 1/2; % smoothness of fluctuations
DEM.M(1).E.d = 2; % embedding dimension
DEM.M(1).E.n = 4; % embedding dimension
DEM.M(1).E.nN = DCM.options.nN; % maximum number of iterations
% adjust M.f (DEM works in time bins not seconds) and initialize M.P
% ---------------------------------------------------------------------
DEM.M(1).delays = M.delays/Y.dt;
DEM.M(1).f = inline([M.f '(x,v,P)*' num2str(Y.dt)],'x','v','P');
% Specify hyper-priors on (log-precision of) observation noise
% ---------------------------------------------------------------------
DEM.M(1).Q = spm_Ce(ones(1,n)); % precision components
DEM.M(1).hE = hE; % prior expectation
DEM.M(1).hC = hC; % prior covariance
% allow (only) neuronal [x, s, f, q, v] hidden states to fluctuate
% ---------------------------------------------------------------------
W = ones(n,1)*[12 16 16 16 16];
DEM.M(1).xP = exp(6); % precision (hidden-state)
DEM.M(1).W = diag(exp(W)); % precision (hidden-motion)
DEM.M(2).V = exp(16); % precision (hidden-cause)
% Generalised filtering (under the Laplace assumption)
% =====================================================================
DEM = spm_LAP(DEM);
% Save DEM estimates
%----------------------------------------------------------------------
DCM.qU = DEM.qU;
DCM.qP = DEM.qP;
DCM.qH = DEM.qH;
% unpack results
% ---------------------------------------------------------------------
F = DEM.F(end);
Ep = DEM.qP.P{1};
Cp = DEM.qP.C;
% predicted responses (y) and residuals (R)
%----------------------------------------------------------------------
y = DEM.qU.v{1}';
R = DEM.qU.z{1}';
R = R - Y.X0*spm_inv(Y.X0'*Y.X0)*(Y.X0'*R);
Ce = exp(-DEM.qH.h{1});
end
% Bilinear representation and first-order hemodynamic kernel
%--------------------------------------------------------------------------
[M0,M1,L1,L2] = spm_bireduce(M,Ep);
[H0,H1] = spm_kernels(M0,M1,L1,L2,M.N,M.dt);
% and neuronal kernels
%--------------------------------------------------------------------------
L = sparse(1:n,(1:n) + 1,1,n,length(M0));
[K0,K1] = spm_kernels(M0,M1,L,M.N,M.dt);
% Bayesian inference and variance {threshold: prior mean plus T = 0}
%--------------------------------------------------------------------------
T = full(spm_vec(pE));
sw = warning('off','SPM:negativeVariance');
Pp = spm_unvec(1 - spm_Ncdf(T,abs(spm_vec(Ep)),diag(Cp)),Ep);
Vp = spm_unvec(full(diag(Cp)),Ep);
warning(sw);
try, M = rmfield(M,'nograph'); end
% Store parameter estimates
%--------------------------------------------------------------------------
DCM.M = M;
DCM.Y = Y;
DCM.U = U;
DCM.Ce = Ce;
DCM.Ep = Ep;
DCM.Cp = Cp;
DCM.Pp = Pp;
DCM.Vp = Vp;
DCM.H1 = H1;
DCM.K1 = K1;
DCM.R = R;
DCM.y = y;
DCM.T = 0;
% Data ID and log-evidence
%==========================================================================
if isfield(M,'FS')
try
ID = spm_data_id(feval(M.FS,Y.y,M));
catch
ID = spm_data_id(feval(M.FS,Y.y));
end
else
ID = spm_data_id(Y.y);
end
% Save approximations to model evidence: negative free energy, AIC, BIC
%--------------------------------------------------------------------------
evidence = spm_dcm_evidence(DCM);
DCM.F = F;
DCM.ID = ID;
DCM.AIC = evidence.aic_overall;
DCM.BIC = evidence.bic_overall;
% Save SPM version and revision number of code used
%--------------------------------------------------------------------------
[DCM.version.SPM.version, DCM.version.SPM.revision] = spm('Ver');
DCM.version.DCM.version = spm_dcm_ui('Version');
DCM.version.DCM.revision = SVNid;
%-Save DCM
%--------------------------------------------------------------------------
if ~isstruct(P)
save(P,'DCM','F','Ep','Cp', spm_get_defaults('mat.format'));
end
if ~nargin
spm('Pointer','Arrow');
spm('FigName','Done');
end
c4science · Help