123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125126127128129130131132133134135136137138139140141142143144145146147148149150151152153154155156157158159160161162163164165166167168169170171172173174175176177178179180181182183184185186187188189190191192193194195196197198199200201202203204205206207208209210211212213214215216217218219220221222223224225226227228229230231232233234235236237238239 |
- function [HMM,csd] = spm_dcm_HMM(GCM,N,b)
- % PEB Inversion of a DCM under a hidden Markov model of state transitions
- % FORMAT [HMM,CSD] = spm_dcm_HMM(DCM,N,b)
- % FORMAT [HMM] = spm_dcm_HMM(CSD,b)
- % -------------------------------------------------------------------------
- % DCM{p} - DCMs for p sessons: DCM.b encodes state-dependent connections
- % N - number of windows within which to evaluate states
- % b{s} - Cell array state transition priors (Dirichlet parameters)
- %
- % returns HMM(s):
- % HMM(s).X - posterior expectation of hidden states
- % HMM(s).qB - posterior expectation of HMM parameters
- % HMM(s).qb - and Dirichlet concentration parameters
- % HMM(s).qP - posterior expectation of PEB parameters
- % HMM(s).qC - posterior covariances of PEB parameters
- % HMM(s).iP - indices of DCM parameters
- % HMM(s).Ep - posterior expectation of DCM parameters
- % HMM(s).Cp - posterior covariances of DCM parameters
- % HMM(s).L - free energy components
- % HMM(s).F - total free energy (model evidence)
- % s - index of HMM structure (prior model of state transitions)
- %
- % CSD{N,P} - inverted DCM of each window; with window functions in CSD{n}.W
- %__________________________________________________________________________
- %
- % This routine characterises a single timeseries in terms of latent or
- % hidden brain states manifest in terms of state dependent connectivity.
- % It first inverts the complex cross spectral density of the observed
- % timeseries and then estimates epoch specific fluctuations in state
- % dependent connectivity in a subset of connections (specified in the
- % logical field DCM.b). The ensuing sequence of posterior densities are
- % then subject to Bayesian model reduction to provide evidence for
- % sequences of state transitions under a hidden Markov model. Effectively,
- % this involves supplying the evidence that the brain is in a particular
- % connectivity state at each epoch – using the reduced free energy - to a
- % variational message passing scheme based upon a Markov decision process.
- % The higher (discrete state space for hidden Markov model) level that
- % returns the Bayesian model average for iterative optimisation of the
- % state dependent connection (PEB) parameters, and the epoch specific
- % connectivity (DCM) parameters. The products of this inversion are
- % posteriors at the DCM (epoch specific), PEB, (state specific)and HMM
- % (transition) level. These posterior densities fully characterise a
- % given time series in terms of discrete state transitions, where each
- % brain state is associated with a location in (connectivity) parameter
- % space; in other words, a discrete characterisation of dynamic or
- % fluctuating effective connectivity.
- %__________________________________________________________________________
- % Copyright (C) 2015-2016 Wellcome Trust Centre for Neuroimaging
- % Karl Friston
- % $Id: spm_dcm_HMM.m 7279 2018-03-10 21:22:44Z karl $
- % get windowed cross spectra if necessary
- %==========================================================================
- if nargin == 3
- csd = spm_dcm_window(GCM,N);
- else
- csd = GCM; % cross spectra of windows
- b = N; % state transition priors
- end
- %% inversion of hierarchical (empirical) Bayesian HMM model
- %==========================================================================
- [N,P] = size(csd); % number of windows and sessions
- a = spm_find_pC(csd{1},{'A'}); % indices of state-dependent parameters
- for k = 1:numel(b)
-
- % initialise this model
- % ---------------------------------------------------------------------
- clear F L M O MDP
- b0 = b{k}; % prior contraints on transitions
- s = size(b0,1); % number of hidden states
- D = sparse(1,1,1,s,1); % intial state
-
- % initialise succession of hidden states (assume an orbit)
- % ---------------------------------------------------------------------
- X = kron(ones(1,N),eye(s,s));
- X = spm_softmax(X(1:s,1:N));
- X = kron(ones(1,P),X);
-
- % initialise priors on hidden Markov model
- % ---------------------------------------------------------------------
- MDP.A = {eye(s,s)}; % likelihood of outcomes given hidden states
- MDP.D = {full(D)}; % initial states
- MDP.b = {full(b0)}; % transitions among states
-
- % prepare
- % ---------------------------------------------------------------------
- for i = 1:16
-
- % update state-dependent parameters; using PEB
- %==================================================================
- M.X = X'; % expected hidden states
- [PEB,CSD] = spm_dcm_peb(csd(:),M,{'A'}); % expected connectivity
-
- % evaluate the likelihood of each state (at each epoch)
- %==================================================================
- for t = 1:numel(CSD)
- qE = spm_vec(CSD{t}.Ep);
- qP = spm_inv(CSD{t}.Cp);
- for o = 1:s
- rE = PEB.Ep(:,o) - qE(a);
- F(o,t) = -rE'*qP(a,a)*rE/2;
- end
- end
- % update expected hidden states, using MDP (HMM)
- %==================================================================
- MDP.O{1} = spm_softmax(F);
- mdp = spm_MDP_VB_X(MDP);
-
- % update transition probabilities and inital states
- %==================================================================
- MDP.B = mdp.b;
- X = mdp.X{1};
-
- % record free energy terms
- % -----------------------------------------------------------------
- L(1,i) = mdp.F(end); % HMM: complexity from hidden states
- L(2,i) = mdp.Fb; % HMM: complexity from HMM parameters
- L(3,i) = PEB.F; % PEB: complexity from PEB parameters
-
-
- % test for convergence
- % -----------------------------------------------------------------
- if i > 4 && norm(L(:,i) - L(:,i - 1)) < 1/128
- break
- end
- end
-
- % store posteriors for this hidden Markov model
- % ---------------------------------------------------------------------
- for t = 1:numel(CSD)
- Ep{t} = CSD{t}.Ep; % epoch specific (PEB) expectation
- Cp{t} = CSD{t}.Cp; % epoch specific (PEB) covariances
- end
-
- HMM(k).X = X; % posterior expectation of hidden states
- HMM(k).qB = mdp.B{1}; % posterior expectation of HMM parameters
- HMM(k).qb = mdp.b{1}; % and Dirichlet concentration parameters
- HMM(k).qP = PEB.Ep; % posterior expectation of PEB parameters
- HMM(k).qC = PEB.Cp; % posterior covariances of PEB parameters
- HMM(k).iP = PEB.Pind; % indices of DCM parameters
- HMM(k).Ep = Ep; % posterior expectation of DCM parameters
- HMM(k).Cp = Cp; % posterior covariances of DCM parameters
- HMM(k).L = L; % free energy components
- HMM(k).F = sum(L(:,end)); % total free energy (model evidence)
-
- end
- if nargout, return, end
- % report analysis
- %==========================================================================
- % get figure
- %--------------------------------------------------------------------------
- spm_figure('Getwin','HMM'); clf
- spm_dcm_HMM_plot(HMM)
- return
- function [csd] = spm_dcm_window(GCM,N)
- % Auxiliary function: cross spectra of windowed segments
- % FORMAT [HMM,CSD] = spm_dcm_HMM(DCM,N)
- % -------------------------------------------------------------------------
- % check DCM is a cell array and number of sessions: e.g., participants (P)
- %--------------------------------------------------------------------------
- if isstruct(GCM), GCM = {GCM}; end
- P = numel(GCM);
- % nonlinear system identification (DCM for CSD)
- %==========================================================================
- % First, invert the entire timeseries to estimate DCM parameters that
- % do not change over time. These parameters will be fixed using precise
- % shrinkage priors to estimate fluctuating (state-dependent) parameters
- % using windowed data to evaluate complex cross spectra
- %--------------------------------------------------------------------------
- % invert sessions and take parameter average
- %--------------------------------------------------------------------------
- for p =1:P
- CSD{p} = spm_dcm_fmri_csd(GCM{p});
- end
- CSD = spm_dcm_bpa(CSD(:));
- % for each session (e.g., subject)
- %--------------------------------------------------------------------------
- for p = 1:P
-
- % preliminaries: extract state-dependent parameter specification (b)
- %----------------------------------------------------------------------
- DCM = GCM{p};
- n = size(DCM.b,1);
- B = DCM.b;
- DCM.b = zeros(n,n,0);
-
- % (overlapping) Hanning windows (epochs) W
- %----------------------------------------------------------------------
- nS = size(DCM.Y.y,1); % number of scans
- nw = nS/N; % window length
- iw = (1:N)*nw - nw/2; % indices of window centres
- jw = -nw:nw;
- hw = spm_hanning(length(jw));
- W = zeros(nS,N);
- for i = 1:N
- k = iw(i) + jw;
- j = find(k > 0 & k < nS);
- W(k(j),i) = hw(j);
- end
-
- % use posterior expectations as precise priors that relax hyperpriors
- %----------------------------------------------------------------------
- pC = spm_zeros(CSD.Ep); % infinitely precise shrinkage priors
- pC.A = B/16; % state dependent parameters
- DCM.M.pE = CSD.Ep; % prior expectation of parameters
- DCM.M.pC = diag(spm_vec(pC)); % prior covariances of parameters
- DCM.M.hE = 0; % expected log degrees of freedom
- DCM.M.hC = 1/64; % intermediate covariance
-
- % invert each window
- %----------------------------------------------------------------------
- DCM.options.order = 6;
- Y.y = DCM.Y.y;
- for t = 1:N
- w = diag(W(:,t)); % get windowing matrix
- w = w(W(:,t) > 0,:); % and truncate timeseries
- DCM.Y.y = w*Y.y; % windowed timeseries
- csd{t,p} = spm_dcm_fmri_csd(DCM); % invert window
- csd{t,p}.W = sum(w)'; % store window
- end
-
- end
|