Homec4science
Diffusion iCAPs public version 2 (master) Browse

iCAPs_v2/packages/spm12/spm_dcm_post_hoc.m
e5f5005bb2efmaster

Last Change
Raw File

spm_dcm_post_hoc.m
View Options

function DCM = spm_dcm_post_hoc(P,fun,field,write_all)
% Post hoc optimisation of DCMs (under the Laplace approximation)
% FORMAT DCM = spm_dcm_post_hoc(P,fun,field,write_all)
%
% P - character/cell array of DCM filenames
% - or cell array of DCM structures; where
% DCM.M.pE - prior expectation (with parameters in pE.A, pE.B and pE.C)
% DCM.M.pC - prior covariance
% DCM.Ep - posterior expectations
% DCM.Cp - posterior covariance
%
% Optional parameters:
%
% fun - optional family definition function: k = fun(A,B,C)
% k = 1,2,...,K for K families or proper subsets of a partition
% of model space - a function of the adjacency matrices: e.g.,
%
% fun = @(A,B,C) any(spm_vec(B(:,:,2))) + 1;
%
% returns 1 if there are no bilinear parameters for the 2nd
% bilinear effect and 2 if there are. fun should be an
% anonymous function or script. NB: Model posteriors over
% families with and without free parameters (in A,B,C and D)
% are evaluated automatically and saved in DCM_BPA (DCM.Pp)
%
% field - the fieldnames of the parameters in the structure pE and Ep
% that are to be included in Bayesian model reduction.
% The default is {'A','B','C'}.
%
% write_all - if true, saves all models from the final iteration of the
% search (i.e. those models in the display) into a subfolder
% named 'reduced' of the original model).
%
%--------------------------------------------------------------------------
% This routine searches over all possible reduced models of a full model
% (DCM) and uses post hoc model selection to select the best. Reduced
% models mean all permutations of free parameters (parameters with a non-
% zero prior covariance), where models are defined in terms of their prior
% covariance. The full model should be inverted prior to post hoc
% optimization. If there are more than 16 free-parameters, this routine
% will implement a greedy search: This entails searching over all
% permutations of the 8 parameters whose removal (shrinking the prior
% variance to zero) produces the smallest reduction (greatest increase)
% in model evidence. This procedure is repeated until all 8 parameters
% are retained in the best model or there are no more parameters to
% consider. When several DCMs are optimized together (as in group studies),
% they are checked to ensure the same free parameters have been specified
% and the log-evidences are pooled in a fixed effects fashion.
%
% This application of post hoc optimization assumes the DCMs that are
% optimized are the same model of different data. Normally, this would be
% a full model, in the sense of having the maximum number of free
% parameters, such that the set of reduced models is as large as possible.
% In contrast spm_dcm_search operates on different DCMs of the same data
% to identify the best model, after inverting the full(est) model
%
% The outputs of this routine are graphics reporting the model reduction
% (optimisation) and a DCM_opt_??? for every specified DCM that contains
% reduced conditional parameters estimates (for simplicity, the original
% kernels and predicted states are retained). The structural and functional
% (spectral embedding) graphs are based on Bayesian parameter averages
% over multiple DCMs, which are stored in DCM_BPA.mat. This DCM also
% contains the posterior probability of models partitioned according to
% whether a particular parameter exists or not:
%
% DCM.Pp - Model posterior (with and without each parameter)
% DCM.Ep - Bayesian parameter average under selected model
% DCM.Cp - Bayesian parameter covariance under selected model
% DCM.Pf - Model posteriors over user specified families
% DCM.fun - User-specified family definition function
% DCM.files - List of DCM files used for Bayesian averaging
%
% See also: spm_dcm_search
%__________________________________________________________________________
% Copyright (C) 2010-2015 Wellcome Trust Centre for Neuroimaging
% Karl Friston, Peter Zeidman
% $Id: spm_dcm_post_hoc.m 6724 2016-02-19 19:13:07Z karl $
%-Number of parameters to consider before invoking greedy search
%--------------------------------------------------------------------------
nmax = 8;
%-Get filenames
%--------------------------------------------------------------------------
if nargin < 1
[P, sts] = spm_select([1 Inf],'^DCM.*\.mat$','Select DCM*.mat files');
if ~sts, return; end
end
if ischar(P), P = cellstr(P); end
if isstruct(P), P = {P}; end
TOL = exp(-16);
%-Check family definition functions
%--------------------------------------------------------------------------
if nargin < 2 || isempty(fun)
fun = [];
end
%-Check fields of parameter stucture
%--------------------------------------------------------------------------
if nargin < 3 || isempty(field)
field = {'A','B','C'};
end
%-Check whether to write all models
%--------------------------------------------------------------------------
if nargin < 4
write_all = false;
else
write_all = logical(write_all);
end
%-Collate these user-supplied parameters
%--------------------------------------------------------------------------
params = struct();
params.fun = fun; % Family function
params.field = field; % Cell array of included fields
params.TOL = TOL; % Numerical tolerance
params.write_all = write_all; % Whether to write all models
params.P = P; % Filenames of all input models
params.N = numel(params.P); % Number of input models
params.nograph = spm('CmdLine'); % Graphical display
%-Run post-hoc DCM
%==========================================================================
% Check that models are compatible in terms of their prior variances
%--------------------------------------------------------------------------
[example_DCM, C] = check_models(params);
% Greedy search (GS) - eliminating parameters in a top down fashion
%--------------------------------------------------------------------------
[C, model_space] = greedy_search(C, example_DCM, nmax, params);
% Inference over families
%--------------------------------------------------------------------------
[Pk, Pf] = family_inference(example_DCM, C, model_space, params);
% Calculate reduced models and BPA
%--------------------------------------------------------------------------
[BPA, P_opt] = compute_post_hoc(example_DCM, C, params);
% Show full and reduced conditional estimates (for Bayesian average)
%--------------------------------------------------------------------------
if ~params.nograph
create_plots(example_DCM, BPA, Pk, Pf, params, ~nargout);
end
% Save Bayesian Parameter Average and family-wise model inference
%--------------------------------------------------------------------------
DCM = save_bpa_dcm(Pk, BPA, Pf, params, P_opt);
%==========================================================================
function [example_DCM,C] = check_models(params)
% Check that models are compatible in terms of their prior variances
%--------------------------------------------------------------------------
% params - user supplied parameters
% example_DCM - an example model
% C - binarized prior variance vector of an example model
for j = 1:params.N
% Get prior covariances
%----------------------------------------------------------------------
try, DCM = load(params.P{j}); DCM = DCM.DCM; catch, DCM = params.P{j}; end
try, pC = diag(DCM.M.pC); catch, pC = spm_vec(DCM.M.pC); end
% and compare it with the first model
%----------------------------------------------------------------------
if j == 1
C = pC;
else
if any(xor(pC,C))
error('Please check model %i for compatibility.\n',j)
end
end
end
example_DCM = DCM;
C = logical(C);
%==========================================================================
function [C,model_space] = greedy_search(C,DCM,nmax,params)
% Perform a greedy search to reduce the number of free parameters in C
%--------------------------------------------------------------------------
% A model space is created, K, each row indicates which free
% parameters to disable. (If there are more than nmax free
% parameters, K will just include combinations of the 8 parameters
% which contribute the least, and the whole process will repeat).
%
% C - binary vector indicating which parameters are free (have
% non-zero prior variance)
% DCM - example model to be reduced
% nmax - number of parameters to consider before invoking greedy
% search
% params - user supplied parameters
%
% C - updated binary vector of all parameters (1=free)
% model_space.k - indices of free parameters to keep
% model_space.K - model space (most recent iteration). 1=disabled parameter
% model_space.p - pooled model evidence (posterior probabilities)
GS = 1;
while GS
%-Find free parameters
%----------------------------------------------------------------------
if isstruct(DCM.Ep)
k = spm_fieldindices(DCM.Ep,params.field{:});
else
k = 1:spm_length(DCM.Ep);
end
k = k(C(k));
nparam = length(k);
%-If there are too many params find those with the least evidence
%----------------------------------------------------------------------
if nparam > nmax
k = identify_parameters_with_least_evidence(k,C,params);
% Flag a greedy search
GS = 1;
elseif isempty(k)
error('There are no free parameters in this model.');
else
GS = 0;
end
%-Evaluate the model space
%----------------------------------------------------------------------
[C,nelim,model_space] = evaluate_model_space(C,k,params);
%-Continue greedy search if any parameters have been eliminated
%----------------------------------------------------------------------
GS = GS & nelim;
%-Show results
%----------------------------------------------------------------------
if ~params.nograph
Fgraph = spm_figure('Getwin','Graphics');
spm_figure('Clear',Fgraph);
fprintf('%i out of %i free parameters removed \n',nelim,nparam)
ax = subplot(3,2,1,'Parent',Fgraph);
if length(model_space.K) > 32, plot(model_space.S,'k'),...
else bar(model_space.S,'c'), end
title(ax,'log-posterior','FontSize',16)
xlabel(ax,'model','FontSize',12)
ylabel(ax,'log-probability','FontSize',12)
axis(ax,'square');
ax = subplot(3,2,2,'Parent',Fgraph);
if length(model_space.K) > 32, plot(model_space.p,'k'),...
else bar(model_space.p,'r'), end
title(ax,'model posterior','FontSize',16)
xlabel(ax,'model','FontSize',12)
ylabel(ax,'probability','FontSize',12)
axis(ax,'square');
drawnow
end
end
%==========================================================================
function [C,nelim,model_space] = evaluate_model_space(C,k,params)
% Evaluate all possible models, each with different permutations of
% parameters (specified in k) eliminated. Return updated C with
% the least successful parameters pruned, plus statistics on the model
% space.
%
% C - binary vector of all parameters (1=free)
% k - indices of free parameters to try disabling
% params - user supplied parameters
%
% C - updated parameter vector (1=free)
% nelim - number of parameters eliminated
% model_space.K - the model space (permutations x parameters in k)
% model_space.k - indices of free parameters (as provided)
% model_space.S - log-posterior probability of each parameter
% model_space.p - posterior probability of each parameter
%-Create model space in terms of free parameter indices
%--------------------------------------------------------------------------
K = spm_perm_mtx(length(k));
%-Loop through DCMs and models and get log-evidences
%--------------------------------------------------------------------------
G = zeros(length(K),params.N);
for j = 1:params.N
try, DCM=load(params.P{j}); DCM=DCM.DCM; catch, DCM = params.P{j}; end
if isstruct(DCM.M.pC), DCM.M.pC = diag(spm_vec(DCM.M.pC)); end
fprintf('\nsearching (%i): 00%%',j)
%-Get priors and posteriors
%----------------------------------------------------------------------
qE = DCM.Ep;
qC = DCM.Cp;
pE = DCM.M.pE;
pC = DCM.M.pC;
%-Remove (a priori) null space
%----------------------------------------------------------------------
U = spm_svd(pC);
qE = U'*spm_vec(qE);
pE = U'*spm_vec(pE);
qC = U'*qC*U;
pC = U'*pC*U;
%-Model search over new prior (covariance)
%----------------------------------------------------------------------
for i = 1:length(K)
r = C; r(k(K(i,:))) = 0;
R = U(r,:)'*U(r,:);
G(i,j) = spm_log_evidence(qE,qC,pE,pC,pE,R*pC*R);
if params.write_all
write_reduced_model(r,i,DCM,params.P{j});
end
fprintf('\b\b\b\b')
fprintf('%-3.0f%%',i*100/length(K))
end
end
fprintf('\n')
%-Pooled model evidence
%--------------------------------------------------------------------------
S = sum(G,2);
S = S - S(end);
p = exp(S - max(S));
p = p/sum(p);
%-Prune parameters (row of K) which gave the largest increase in
% evidence when removed
% -------------------------------------------------------------------------
[unused, i] = max(p);
C(k(K(i,:))) = 0;
% Number of eliminated parameters
%--------------------------------------------------------------------------
nelim = full(sum(K(i,:)));
model_space = struct('K',K, 'k',k, 'S',S, 'p',p);
%==========================================================================
function k = identify_parameters_with_least_evidence(k,C,params)
% Identify the 8 parameters with the least evidence
%--------------------------------------------------------------------------
% k - Updated indices of free parameters (with least evidence)
% C - Binary vector of all parameters (1=free)
% params - User supplied parameters
%
% k - Indices of free parameters
%-Loop through DCMs and free parameters and get log-evidences
%--------------------------------------------------------------------------
Z = zeros(length(k),params.N);
for j = 1:params.N
try, DCM=load(params.P{j}); DCM=DCM.DCM; catch, DCM = params.P{j}; end
if isstruct(DCM.M.pC), DCM.M.pC = diag(spm_vec(DCM.M.pC)); end
fprintf('\ninitial search (%i): 00%%',j)
%-Get priors and posteriors
%----------------------------------------------------------------------
qE = DCM.Ep;
qC = DCM.Cp;
pE = DCM.M.pE;
pC = DCM.M.pC;
%-Remove (a priori) null space
%----------------------------------------------------------------------
U = spm_svd(pC);
qE = U'*spm_vec(qE);
pE = U'*spm_vec(pE);
qC = U'*qC*U;
pC = U'*pC*U;
%-Model search over new prior without the i-th parameter
%----------------------------------------------------------------------
for i = 1:length(k)
r = C; r(k(i)) = 0;
R = U(r,:)'*U(r,:);
Z(i,j) = spm_log_evidence(qE,qC,pE,pC,pE,R*pC*R);
fprintf('\b\b\b\b')
fprintf('%-3.0f%%',i*100/length(k))
end
end
%-Find parameters with the least evidence
%--------------------------------------------------------------------------
Z = sum(Z,2);
[unused,i] = sort(-Z);
k = k(i(1:8));
%==========================================================================
function [Pk, Pf] = family_inference(DCM,C,model_space,params)
% Perform family level inference
%--------------------------------------------------------------------------
% DCM - example DCM structure
% C - binary vector of all free parameters
% model_space.K - model space
% model_space.k - vector of free parameter indices
% model_space.p - posterior probabiliy of all models
% params - User supplied parameters
%
% Pk - posterior probability of each parameter
% Pf - posterior probability of each family
pE = DCM.M.pE;
%-Inference over families (one family per coupling parameter)
%--------------------------------------------------------------------------
Pk = zeros(2,length(model_space.k));
for i = 1:length(model_space.k)
Pk(1,i) = mean(model_space.p(~model_space.K(:,i)));
Pk(2,i) = mean(model_space.p( model_space.K(:,i)));
end
Pk = Pk(1,:)./sum(Pk);
Pn = full(double(C));
Pn(model_space.k) = Pk;
Pk = spm_unvec(Pn,pE);
%-Inference over families (specified by fun)
%--------------------------------------------------------------------------
if isempty(params.fun)
Pf = 0;
else
for i = 1:length(model_space.K)
Pn = full(double(C));
Pn(model_space.K(i,:)) = 0;
Pn = spm_unvec(Pn,pE);
try
Kf(i) = params.fun(Pn.A,Pn.B,Pn.C);
catch
try
Kf(i) = params.fun(Pn.A,Pn.B);
catch
Kf(i) = params.fun(Pn.A);
end
end
end
for i = 1:max(Kf)
Pf(i) = mean(model_space.p(Kf == i));
end
Pf(isnan(Pf)) = 0;
Pf = Pf/sum(Pf);
end
%==========================================================================
function [BPA,P_opt] = compute_post_hoc(DCM,C,params)
% For each model, create reduced model (DCM_opt*.mat) and collate data
% for bayesian parameter average (CQ,Cq,EQ,Eq).
%--------------------------------------------------------------------------
% DCM - template DCM structure
% C - binary vector of parameters
% params - User supplied parameters
%
% BPA.CQ - BPA covariances (full model)
% BPA.Cq - BPA covariances (reduced model)
% BPA.EQ - BPA means (full model)
% BPA.Eq - BPA means (reduced model)
% P_opt - Filenames of optimal models
if isstruct(DCM.M.pC), DCM.M.pC = diag(spm_vec(DCM.M.pC)); end
%-Selected model (reduced prior covariance)
%--------------------------------------------------------------------------
pC = DCM.M.pC;
rC = diag(C)*pC*diag(C);
R = spm_unvec(full(C),DCM.M.pE);
try
R.a = DCM.a & R.A;
R.b = DCM.b & R.B;
R.c = DCM.c & R.C;
R.d = DCM.d & R.D;
end
EQ = 0;
PQ = 0;
Eq = 0;
Pq = 0;
P_opt = {};
for j = 1:params.N
%-Get priors and posteriors
%----------------------------------------------------------------------
try, DCM=load(params.P{j}); DCM=DCM.DCM; catch, DCM = params.P{j}; end
qE = DCM.Ep;
qC = DCM.Cp;
pE = DCM.M.pE;
pC = DCM.M.pC;
%-Get posterior of selected model - rC
%----------------------------------------------------------------------
[F,Ep,Cp] = spm_log_evidence_reduce(qE,qC,pE,pC,pE,rC);
%-Bayesian parameter average (for full and reduced selected model)
%----------------------------------------------------------------------
PP = spm_inv(qC,params.TOL);
Pp = spm_inv(Cp,params.TOL);
EQ = EQ + PP*spm_vec(qE); % BPA means (full model)
Eq = Eq + Pp*spm_vec(Ep); % BPA means (reduced model)
PQ = PQ + PP; % BPA precisions (full model)
Pq = Pq + Pp; % BPA precisions (reduced model)
%-Put reduced conditional estimates in DCM
%----------------------------------------------------------------------
% Bayesian inference and variance
try
T = full(spm_vec(pE)) + DCM.T;
catch
T = full(spm_vec(pE));
end
sw = warning('off','SPM:negativeVariance');
Pp = spm_unvec(1 - spm_Ncdf(T,abs(spm_vec(Ep)),diag(Cp)),Ep);
Vp = spm_unvec(diag(Cp),Ep);
warning(sw);
% Store parameter estimates
DCM.M.pC = rC;
DCM.Ep = Ep;
DCM.Cp = Cp;
DCM.Pp = Pp;
DCM.Vp = Vp;
% and in DEM format
DCM.qP.P{1} = Ep;
DCM.qP.C = Cp;
DCM.qP.V{1} = spm_unvec(diag(Cp),Ep);
% and in prior constraints fields
try
DCM.a = R.a;
DCM.b = R.b;
DCM.c = R.c;
DCM.d = R.d;
end
% approximations to model evidence: negative free energy, AIC, BIC
DCM.F = F;
try
evidence = spm_dcm_evidence(DCM);
DCM.AIC = evidence.aic_overall;
DCM.BIC = evidence.bic_overall;
end
%-Save optimised DCM
%------------------------------------------------------------------
try
[pth, name] = fileparts(params.P{j});
if ~strncmp(name,'DCM_opt_',8)
name = ['DCM_opt_' name(5:end) '.mat'];
end
P_opt{j} = fullfile(pth,name);
catch
P_opt{j} = fullfile(pwd,['DCM_opt_' date '.mat']);
end
save(P_opt{j},'DCM','F','Ep','Cp', spm_get_defaults('mat.format'));
end
%-Calculate bayesian parameter average
%--------------------------------------------------------------------------
if isstruct(pC), pC = diag(spm_vec(pC)); end
CQ = spm_inv(PQ,params.TOL); % post covariance (full)
Cq = spm_inv(Pq,params.TOL); % post covariance (reduced)
EQ = CQ*EQ; % post mean (full)
Eq = Cq*Eq; % post mean (reduced)
EQ = spm_unvec(EQ,pE);
Eq = spm_unvec(Eq,pE);
% Pack BPA results to return
%--------------------------------------------------------------------------
BPA.CQ = CQ; BPA.Cq = Cq;
BPA.EQ = EQ; BPA.Eq = Eq;
%==========================================================================
function create_plots(DCM,BPA,Pk,Pf,params,extra_plots)
% Plot results
%--------------------------------------------------------------------------
% DCM - Template DCM structure
% BPA.CQ - BPA posterior covariance (full)
% BPA.Cq - BPA posterior covariance (reduced)
% BPA.EQ - BPA posterior mean (full)
% BPA.Eq - BPA posterior mean (reduced)
% Pk - Posterior probability of each parameter
% Pf - Posterior probability of each family
% params - User supplied parameters
% extra_plots - additional plots if there are no output arguments
if isstruct(DCM.Ep)
i = spm_fieldindices(DCM.Ep,'A','B','C');
else
i = 1:spm_length(DCM.Ep);
end
pE = DCM.M.pE;
pE = spm_vec(pE);
EP = spm_vec(BPA.EQ);
Ep = spm_vec(BPA.Eq);
%-Plot summary of results
%--------------------------------------------------------------------------
Fgraph = spm_figure('GetWin','Graphics');
ax = subplot(3,2,3,'Parent',Fgraph);
spm_plot_ci(EP(i),BPA.CQ(i,i));
title(ax,'MAP connections (full)','FontSize',16);
axis(ax,'square');
ax = subplot(3,2,4,'Parent',Fgraph);
spm_plot_ci(Ep(i),abs(BPA.Cq(i,i)));
title(ax,'MAP connections (reduced)','FontSize',16);
axis(ax,'square');
ax = subplot(3,2,5,'Parent',Fgraph);
bar(ax,EP(i) - pE(i));
xlabel(ax,'parameter');
title(ax,'MAP minus prior','FontSize',16);
spm_axis(ax,'tight');
axis(ax,'square');
ax = subplot(3,2,6,'Parent',Fgraph);
bar(ax,Ep(i) - EP(i));
xlabel(ax,'parameter');
title(ax,'differences in MAP','FontSize',16);
spm_axis(ax,'tight');
axis(ax,'square');
drawnow
%-Stop this function here if we only need basic plots
%--------------------------------------------------------------------------
if ~extra_plots
return;
end
%-Show structural and functional graphs
%--------------------------------------------------------------------------
F = spm_figure('GetWin','Graph analysis');
spm_figure('Clear',F);
try
spm_dcm_graph(DCM.xY,BPA.Eq.A);
catch
try
spm_dcm_graph(DCM,BPA.Eq.A);
catch
delete(F);
end
end
%-Show coupling matrices
%--------------------------------------------------------------------------
if numel(params.field) == 3 && isfield(DCM,'a')
F = spm_figure('GetWin','Bayesian parameter average (selected model)');
spm_figure('Clear',F);
spm_dcm_fmri_image(BPA.Eq);
F = spm_figure('GetWin','Model posterior (over parameters)');
spm_figure('Clear',F);
spm_dcm_fmri_image(Pk);
end
%-Illustrate family-wise inference
%--------------------------------------------------------------------------
if ~isempty(params.fun)
F = spm_figure('GetWin','Model posterior (over families)');
ax = subplot(2,1,1,'Parent',F);
bar(ax,Pf);
xlabel(ax,'family');
title(ax,'Model posterior (over families)','FontSize',16);
axis(ax,'square');
end
%==========================================================================
function DCM = save_bpa_dcm(Pk,BPA,Pf,params,P_opt)
% Save the Bayesian Parameter Average
%--------------------------------------------------------------------------
% Pk - Posterior probability of each parameter
% BPA.Cq - BPA covariances (reduced model)
% BPA.Eq - BPA means (reduced model)
% Pf - Model posteriors over user specified families
% params - User supplied parameters
% P_opt - Optimal model filenames
% Get original (first) DCM
%--------------------------------------------------------------------------
try DCM=load(P_opt{1}); DCM=DCM.DCM; catch, DCM = P_opt{1}; end
DCM.Pp = Pk;
DCM.Ep = BPA.Eq;
DCM.Cp = BPA.Cq;
DCM.fun = params.fun;
DCM.files = P_opt;
DCM.Pf = Pf;
% and save as DCM_BPA
%--------------------------------------------------------------------------
try
pth = fileparts(P_opt{1});
name = 'DCM_BPA.mat';
name = fullfile(pth,name);
catch
name = fullfile(pwd,'DCM_BPA.mat');
end
save(name,'DCM', spm_get_defaults('mat.format'));
%==========================================================================
function write_reduced_model(C,i,DCM,filename)
% Create and save a DCM file for a reduced model
%--------------------------------------------------------------------------
% C - reduced variance vector
% i - model index
% DCM - model
% filename - original filename of the model
R = spm_unvec(full(C),DCM.M.pE);
try
R.a = DCM.a & R.A;
R.b = DCM.b & R.B;
R.c = DCM.c & R.C;
R.d = DCM.d & R.D;
end
%-Get priors and posteriors
%--------------------------------------------------------------------------
qE = DCM.Ep;
qC = DCM.Cp;
pE = DCM.M.pE;
pC = DCM.M.pC;
rC = diag(C)*pC*diag(C);
%-Get posterior of selected model - rC
%--------------------------------------------------------------------------
[F,Ep,Cp] = spm_log_evidence_reduce(qE,qC,pE,pC,pE,rC);
%-Put reduced conditional estimates in DCM
%--------------------------------------------------------------------------
% Bayesian inference and variance
%--------------------------------------------------------------------------
try
T = full(spm_vec(pE)) + DCM.T;
catch
T = full(spm_vec(pE));
end
sw = warning('off','SPM:negativeVariance');
Pp = spm_unvec(1 - spm_Ncdf(T,abs(spm_vec(Ep)),diag(Cp)),Ep);
Vp = spm_unvec(diag(Cp),Ep);
warning(sw);
% Store parameter estimates
%--------------------------------------------------------------------------
DCM.M.pC = rC;
DCM.Ep = Ep;
DCM.Cp = Cp;
DCM.Pp = Pp;
DCM.Vp = Vp;
% and in DEM format
%--------------------------------------------------------------------------
DCM.qP.P{1} = Ep;
DCM.qP.C = Cp;
DCM.qP.V{1} = spm_unvec(diag(Cp),Ep);
% and in prior constraints fields
%--------------------------------------------------------------------------
try
DCM.a = R.a;
DCM.b = R.b;
DCM.c = R.c;
DCM.d = R.d;
end
% approximations to model evidence: negative free energy, AIC, BIC
%--------------------------------------------------------------------------
DCM.F = F;
try
evidence = spm_dcm_evidence(DCM);
DCM.AIC = evidence.aic_overall;
DCM.BIC = evidence.bic_overall;
end
%-Save as DCM_reduced_*.mat
%--------------------------------------------------------------------------
try
[pth,name] = fileparts(filename);
pth = fullfile(pth,'reduced');
if ~exist(pth,'file')
mkdir(pth);
end
name = sprintf('%s_reduced_%04i.mat',name,i);
name = fullfile(pth,name);
catch
name = sprintf('DCM_reduced_%04i.mat',i);
end
save(name,'DCM', spm_get_defaults('mat.format'));
c4science · Help