Homec4science
Diffusion iCAPs public (master) Browse

iCAPs/packages/spm12/spm_ar_reml.m
31022b3a13e6master

Last Change
Raw File

spm_ar_reml.m
View Options

function [C,h,Ph,F] = spm_ar_reml(YY,X,m,N)
% ReML estimation of covariance components from y*y'
% FORMAT [C,h,Ph,F] = spm_ar_reml(YY,X,m,N);
%
% YY - (m x m) sample covariance matrix Y*Y' {Y = (m x N) data matrix}
% X - (m x p) design matrix
% m - (1) order of AR(m) model
% N - number of samples
%
%
% C - (m x m) estimated errors = h(1)*Q{1} + h(2)*Q{2} + ...
% h - (q x 1) ReML hyperparameters h: normalised AR coeficients
% Ph - (q x q) conditional precision of h (unnormalised)
%
% F - [-ve] free energy F = log evidence = p(Y|X,Q) = ReML objective
%
% Performs a Fisher-Scoring ascent on F to find ReML variance parameter
% estimates.
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% John Ashburner & Karl Friston
% $Id: spm_ar_reml.m 5219 2013-01-29 17:07:07Z spm $
% assume a single sample if not specified
%--------------------------------------------------------------------------
try
N;
catch
N = 1;
end
% assume AR(1) if not specified
%--------------------------------------------------------------------------
try
m;
catch
m = 1;
end
% ortho-normalise X
%--------------------------------------------------------------------------
if isempty(X)
X = sparse(length(YY),1);
else
X = orth(full(X));
end
% initialise h
%--------------------------------------------------------------------------
m = m + 1;
n = length(YY);
dh = zeros(m,1);
dFdh = zeros(m,1);
dFdhh = zeros(m,m);
L = zeros(m,m);
% initialise and specify hyperpriors
%--------------------------------------------------------------------------
hE = sparse(1,1,1,m,1);
hP = speye(m,m)/exp(32);
h = hE;
% initialise precision components
%--------------------------------------------------------------------------
for i = 1:m
Q{i} = spdiags(ones(n,2),[(1 - i) (i - 1)],n,n);
end
% scale data
%--------------------------------------------------------------------------
Ys = norm(YY,1)/N;
YY = YY/Ys;
% ReML (EM/VB)
%--------------------------------------------------------------------------
for k = 1:64
% compute current estimate of covariance
%----------------------------------------------------------------------
iC = sparse(n,n);
for i = 1:m
iC = iC + Q{i}*h(i);
end
C = inv(iC);
% E-step: conditional covariance cov(B|y) {Cq}
%======================================================================
Cq = pinv(X'*iC*X);
% M-step: ReML estimate of hyperparameters
%======================================================================
% Gradient dF/dh (first derivatives)
%----------------------------------------------------------------------
P = C - X*Cq*X';
U = (iC*YY/N*iC - iC)*P;
for i = 1:m
% dF/dh = -trace(dF/diC*iC*Q{i}*iC)
%------------------------------------------------------------------
QP{i} = Q{i}*P;
dFdh(i) = -trace(QP{i}*U)*N/2;
end
% Expected curvature E{dF/dhh} (second derivatives)
%----------------------------------------------------------------------
for i = 1:m
for j = i:m
% dF/dhh = -trace{P*Q{i}*P*Q{j}}
%--------------------------------------------------------------
dFdhh(i,j) = -trace(QP{i}*QP{j})*N/2;
dFdhh(j,i) = dFdhh(i,j);
end
end
% add hyperpriors
%----------------------------------------------------------------------
e = h - hE;
dFdh = dFdh - hP*e;
dFdhh = dFdhh - hP;
% update regulariser
%----------------------------------------------------------------------
if ~rem(k,8)
L = speye(m,m)*norm(dFdhh)/128;
end
% Fisher scoring: update dh = -inv(ddF/dhh)*dF/dh
%----------------------------------------------------------------------
Ph = -dFdhh;
dh = -pinv(dFdhh - L)*dFdh;
% preclude numerical overflow
%----------------------------------------------------------------------
h = h + dh;
% Convergence (1% change in log-evidence)
%======================================================================
dF = dFdh'*dh;
fprintf('%-30s: %i %30s%e\n',' ReML Iteration',k,'...',full(dF));
if dF < 1e-1, break, end
end
% log evidence = ln p(y|X,Q) = ReML objective = F = trace(R'*iC*R*YY)/2 ...
%--------------------------------------------------------------------------
if nargout > 3
R = P*iC;
F = - trace(R*YY*R')/2 ...
- e'*hP*e/2 ...
- N*n*log(2*pi)/2 ...
- N*spm_logdet(C)/2 ...
+ N*spm_logdet(Cq)/2 ...
- spm_logdet(Ph)/2 ...
+ spm_logdet(hP)/2;
end
% rescale (NB - Q{1) = 2*I
%--------------------------------------------------------------------------
C = C*Ys;
h = -h(2:m)/(h(1)*2);
c4science · Help