123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657585960616263646566676869707172737475767778798081828384858687888990919293949596979899100101102103104105106107108109110111112113114115116117118119120121122123124125 |
- function model = spm_mvb_G(X,L,X0,G,V)
- % Multivariate Bayesian inversion of a linear model
- % FORMAT model = spm_mvb_G(X,L,X0,G,V);
- % X - contrast or target vector
- % L - pattern matrix (n x m)
- % X0 - confounds
- % G - pattern subsets (in columns of G) (m x h)
- % V - cell array of observation noise covariance components
- %
- % returns model:
- % F: log-evidence [F(0), F(1),...]
- % G: pattern switches
- % h: covariance hyperparameters (on R and cov(E))
- % qE: conditional expectation of pattern-weights
- % MAP: MAP projector (pattern-weights)
- % Cp: prior covariance (pattern space)
- %__________________________________________________________________________
- %
- % model: X = L*P + X0*Q + R
- % P = E;
- % cov(E) = h1*diag(G(:,1)) + h2*diag(G(:,2)) + ...
- %
- % See spm_mvb and:
- %
- % Bayesian decoding of brain images.
- % Friston K, Chu C, Mourão-Miranda J, Hulme O, Rees G, Penny W, Ashburner J.
- % Neuroimage. 2008 Jan 1;39(1):181-205
- %
- % Multiple sparse priors for the M/EEG inverse problem.
- % Friston K, Harrison L, Daunizeau J, Kiebel S, Phillips C, Trujillo-Barreto
- % N, Henson R, Flandin G, Mattout J.
- % Neuroimage. 2008 Feb 1;39(3):1104-20.
- %
- % Characterizing dynamic brain responses with fMRI: a multivariate approach.
- % Friston KJ, Frith CD, Frackowiak RS, Turner R.
- % Neuroimage. 1995 Jun;2(2):166-72.
- %__________________________________________________________________________
- % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
-
- % Karl Friston
- % $Id: spm_mvb_G.m 3813 2010-04-07 19:21:49Z karl $
-
- % defaults
- %--------------------------------------------------------------------------
- Nx = size(X,1);
- Nk = size(G,1);
- Np = size(G,2);
- if isempty(X0), X0 = zeros(Nx,1); end
- try, V; catch, V = speye(Nx); end
-
- % null space of confounds
- %--------------------------------------------------------------------------
- X0 = full(X0);
- R = spm_svd(speye(size(X0,1)) - X0*pinv(X0));
- L = R'*L;
- X = R'*X;
- % random effects (and serial correlations)
- %--------------------------------------------------------------------------
- if iscell(V)
- Ne = length(V);
- for i = 1:Ne
- Qe{i} = R'*V{i}*R;
- end
- elseif isnumeric(V)
- Ne = 1;
- Qe = {R'*V*R};
- else
- Ne = 1;
- Qe = {speye(R'*R)};
- end
-
- % assemble empirical priors
- %==========================================================================
- Qp = {};
- LQpL = {};
- for i = 1:Np
- j = find(G(:,i));
- Qp{i} = sparse(j,j,1,Nk,Nk);
- LQpL{i} = L*Qp{i}*L';
- end
-
- % inverse solution
- %==========================================================================
-
- % Covariance components (with the first error Qe{1} fixed)
- %--------------------------------------------------------------------------
- if size(L,2) > 0
- Q = [Qe LQpL];
- else
- Q = Qe;
- end
- % ReML (with mildly informative hyperpriors) and a lower bound on error
- %--------------------------------------------------------------------------
- N = size(X,2);
- YY = X*X';
- Qe = exp(-2)*trace(YY)*Q{1}/trace(Q{1})/N;
- hE = -4;
- hC = 16;
- [Cy,h,P,F] = spm_reml_sc(YY,[],Q,N,hE,hC,Qe);
- % prior covariance: source space
- %--------------------------------------------------------------------------
- Cp = sparse(Nk,Nk);
- for i = 1:Np
- Cp = Cp + h(i + Ne)*Qp{i};
- end
-
- % MAP estimates of instantaneous sources
- %==========================================================================
- MAP = Cp*L'/Cy*R';
- qE = MAP*R*X;
-
-
- % assemble results (pattern weights)
- %==========================================================================
- model.F = F;
- model.G = G;
- model.h = h;
- model.qE = qE;
- model.MAP = MAP;
- model.Cp = Cp;
|