In [35]:
# TO DO WITH THIS FUNC:
# Try it working up to spm_diff
# ...checking vs. spm func using spm debugger breakpoints
# then come up with alternative for spm_diff
# then can do full simulations
In [187]:
import scipy.sparse as sps
code
In [405]:
def fx_erp(x,u,P,M):
"""
% State equations for a neural mass model of erps
Port of spm_fx_erp.m
Usage:
f,J,D = spm_fx_erp(x,u,P,M,returnJ=True,returnD=True)
x - state vector
x(:,1) - voltage (spiny stellate cells)
x(:,2) - voltage (pyramidal cells) +ve
x(:,3) - voltage (pyramidal cells) -ve
x(:,4) - current (spiny stellate cells) depolarizing
x(:,5) - current (pyramidal cells) depolarizing
x(:,6) - current (pyramidal cells) hyperpolarizing
x(:,7) - voltage (inhibitory interneurons)
x(:,8) - current (inhibitory interneurons) depolarizing
x(:,9) - voltage (pyramidal cells)
f - dx(t)/dt = f(x(t))
J - df(t)/dx(t)
D - delay operator dx(t)/dt = f(x(t - d)) = D(d)*f(x(t))
Prior fixed parameter scaling [Defaults]
M.pF.E = [32 16 4]; % extrinsic rates (forward, backward, lateral)
M.pF.H = [1 4/5 1/4 1/4]*128; % intrinsic rates (g1, g2 g3, g4)
M.pF.D = [2 16]; % propogation delays (intrinsic, extrinsic)
M.pF.G = [4 32]; % receptor densities (excitatory, inhibitory)
M.pF.T = [8 16]; % synaptic constants (excitatory, inhibitory)
M.pF.S = [1 1/2]; % parameters of activation function
JG 2017
"""
"""
% get dimensions and configure state variables
%--------------------------------------------------------------------------
n = length(P.A{1}); % number of sources
x = spm_unvec(x,M.x); % neuronal states
""";
n = len(P['A'][0]); # % number of sources # TO DO
x = spm_unvec(x,M['x']); # % neuronal states
"""
#% [default] fixed parameters
#%--------------------------------------------------------------------------
E = [1 1/2 1/8]*32; % extrinsic rates (forward, backward, lateral)
G = [1 4/5 1/4 1/4]*128; % intrinsic rates (g1 g2 g3 g4)
D = [2 16]; % propogation delays (intrinsic, extrinsic)
H = [4 32]; % receptor densities (excitatory, inhibitory)
T = [8 16]; % synaptic constants (excitatory, inhibitory)
R = [2 1]/3; % parameters of static nonlinearity
""";
E = np.array([1., 1/2., 1/8.,])*32; # extrinsic rates (forward, backward, lateral)
G = np.array([1, 4/5., 1/4., 1/4*128.]); # % intrinsic rates (g1 g2 g3 g4)
D = np.array([2, 16]); # % propogation delays (intrinsic, extrinsic)
H = np.array([4, 32]); # % receptor densities (excitatory, inhibitory)
T = np.array([8, 16]); # % synaptic constants (excitatory, inhibitory)
R = np.array([2, 1])/3.; #% parameters of static nonlinearity
"""
% [specified] fixed parameters
%--------------------------------------------------------------------------
if isfield(M,'pF')
try, E = M.pF.E; end
try, G = M.pF.H; end
try, D = M.pF.D; end
try, H = M.pF.G; end
try, T = M.pF.T; end
try, R = M.pF.R; end
end
"""
if 'pF' in M: # if isfield(M,'pF')
try: E = M['pF']['E'];
except: _
try: G = M['pF']['H'];
except: _
try: D = M['pF']['D'];
except: _
try: H = M['pF']['G'];
except: _
try: T = M['pF']['T'];
except: _
try: R = M['pF']['R'];
except: _
"""
#% test for free parameters on intrinsic connections
#%--------------------------------------------------------------------------
try
G = G.*exp(P.H);
end
G = ones(n,1)*G;
"""
try:
G = G*np.exp(P['H']);
except: _
G = np.ones([n,1])*G;
"""
% exponential transform to ensure positivity constraints
%--------------------------------------------------------------------------
A{1} = exp(P.A{1})*E(1);
A{2} = exp(P.A{2})*E(2);
A{3} = exp(P.A{3})*E(3);
C = exp(P.C);
"""
A[0] = np.exp(P['A'][0])*E[0]
A[1] = np.exp(P['A'][1])*E[1]
A[2] = np.exp(P['A'][2])*E[2]
C = np.exp(P['C'])
"""
% intrinsic connectivity and parameters
%--------------------------------------------------------------------------
Te = T(1)/1000*exp(P.T(:,1)); % excitatory time constants
Ti = T(2)/1000*exp(P.T(:,2)); % inhibitory time constants
He = H(1)*exp(P.G(:,1)); % excitatory receptor density
Hi = H(2)*exp(P.G(:,2)); % inhibitory receptor density
"""
Te = T[0]/1000*np.exp(P['T'][:,0])
Ti = T[1]/1000*np.exp(P['T'][:,1])
He = H[0]*np.exp(P['G'][:,0])
Hi = H[1]*np.exp(P['G'][:,1])
"""
% pre-synaptic inputs: s(V)
%--------------------------------------------------------------------------
R = R.*exp(P.S);
S = 1./(1 + exp(-R(1)*(x - R(2)))) - 1./(1 + exp(R(1)*R(2)));
"""
R = R*np.exp(P['S'])
S = 1./(1 + np.exp(-R[0]*(x - R[1]))) - 1./(1 + np.exp(R[0]*R[1]));
"""
% input
%==========================================================================
if isfield(M,'u')
% endogenous input
%----------------------------------------------------------------------
U = u(:)*64;
else
% exogenous input
%----------------------------------------------------------------------
U = C*u(:)*2;
end
"""
if 'u' in M:
#% endogenous input
#%----------------------------------------------------------------------
#U = u(:)*64;
U = u[:]*64.
else:
#% exogenous input
#%----------------------------------------------------------------------
#U = C*u(:)*2;
U = C*u[:]*2
#% State: f(x)
#%==========================================================================
"""
% Supragranular layer (inhibitory interneurons): Voltage & depolarizing current
%--------------------------------------------------------------------------
f(:,7) = x(:,8);
f(:,8) = (He.*((A{2} + A{3})*S(:,9) + G(:,3).*S(:,9)) - 2*x(:,8) - x(:,7)./Te)./Te;
"""
f[:,6] = x[:,7]
f[:,7] = (He*( (A[1] + A[2]) * S[:,8] + G[:,2]*S[:,8]) - 2*x[:,7] - x[:,6]/Te)/Te
"""
% Granular layer (spiny stellate cells): Voltage & depolarizing current
%--------------------------------------------------------------------------
f(:,1) = x(:,4);
f(:,4) = (He.*((A{1} + A{3})*S(:,9) + G(:,1).*S(:,9) + U) - 2*x(:,4) - x(:,1)./Te)./Te;
"""
f[:,0] = x[:,3];
f[:,3] = (He*((A[0] + A[2])*S[:,8] + G[:,0]*S[:,8] + U) - 2*x[:,3] - x[:,0]/Te)/Te;
"""
% Infra-granular layer (pyramidal cells): depolarizing current
%--------------------------------------------------------------------------
f(:,2) = x(:,5);
f(:,5) = (He.*((A{2} + A{3})*S(:,9) + G(:,2).*S(:,1)) - 2*x(:,5) - x(:,2)./Te)./Te;
"""
f[:,1] = x[:,4];
f[:,4] = (He*((A[1] + A[2])*S[:,8] + G[:,1]*S[:,0]) - 2*x[:,4] - x[:,1]/Te)/Te;
"""
% Infra-granular layer (pyramidal cells): hyperpolarizing current
%--------------------------------------------------------------------------
f(:,3) = x(:,6);
f(:,6) = (Hi.*G(:,4).*S(:,7) - 2*x(:,6) - x(:,3)./Ti)./Ti;
"""
f[:,2] = x[:,5];
f[:,5] = (Hi*G[:,3]*S[:,6] - 2*x[:,5] - x[:,2]/Ti)/Ti;
"""
% Infra-granular layer (pyramidal cells): Voltage
%--------------------------------------------------------------------------
f(:,9) = x(:,5) - x(:,6);
f = spm_vec(f);
"""
f[:,8] = x[:,4] - x[:,5];
#f = spm_vec(f);
f = spm_vec(f)
#if nargout < 2; return, end
if returnJ == False and returnD == False:
return f
else:
"""
% Jacobian
%==========================================================================
J = spm_diff(M.f,x,u,P,M,1);
"""
# TO DO
#J = spm_diff(M.f,x,u,P,M,1);
J = diff_dfdx(M['f'],x,u,P,M,1)
"""
% delays
%==========================================================================
% Delay differential equations can be integrated efficiently (but
% approximately) by absorbing the delay operator into the Jacobian
%
% dx(t)/dt = f(x(t - d))
% = Q(d)f(x(t))
%
% J(d) = Q(d)df/dx
%--------------------------------------------------------------------------
De = D(2).*exp(P.D)/1000;
Di = D(1)/1000;
De = (1 - speye(n,n)).*De;
Di = (1 - speye(9,9)).*Di;
De = kron(ones(9,9),De);
Di = kron(Di,speye(n,n));
D = Di + De;
"""
De = D[1] * np.exp(P['D'])/1000.
Di = D[0] / 1000.
De = np.array*([1-sps.eye(n,n)])*De
Di = np.array*([1-sps.eye(9,9)])*Di
De = np.kron(np.ones([9,9]),De)
Di = np.kron(Di,sps.eye([n,n]))
D = Di + De
"""
% Implement: dx(t)/dt = f(x(t - d)) = inv(1 + D.*dfdx)*f(x(t))
% = Q*f = Q*J*x(t)
%--------------------------------------------------------------------------
Q = inv(speye(length(J)) + D.*J);
""";
Q = np.linalg.inv(sps.eye(len(J)) + D*J)
return f,j,D
tests:
In [409]:
%load /home/jgriffiths/Code/libraries_of_others/misc/SPM/spm12/spm_dfdx.m
In [ ]:
# %load /home/jgriffiths/Code/libraries_of_others/misc/SPM/spm12/spm_diff.m
In [ ]:
def spm_diff(f,x,n,V=None,q=None):
"""
matrix high-order numerical differnetiation
port of spm_diff.m
JG June 2017
"""
% matrix high-order numerical differentiation
% FORMAT [dfdx] = spm_diff(f,x,...,n)
% FORMAT [dfdx] = spm_diff(f,x,...,n,V)
% FORMAT [dfdx] = spm_diff(f,x,...,n,'q')
%
% f - [inline] function f(x{1},...)
% x - input argument[s]
% n - arguments to differentiate w.r.t.
%
% V - cell array of matrices that allow for differentiation w.r.t.
% to a linear transformation of the parameters: i.e., returns
%
% df/dy{i}; x = V{i}y{i}; V = dx(i)/dy(i)
%
% q - (char) flag to preclude default concatenation of dfdx
%
% dfdx - df/dx{i} ; n = i
% dfdx{p}...{q} - df/dx{i}dx{j}(q)...dx{k}(p) ; n = [i j ... k]
%
%
% - a cunning recursive routine
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_diff.m 6110 2014-07-21 09:36:13Z karl $
% create inline object
%--------------------------------------------------------------------------
f = spm_funcheck(varargin{1});
% parse input arguments
%--------------------------------------------------------------------------
if iscell(varargin{end})
x = varargin(2:(end - 2));
n = varargin{end - 1};
V = varargin{end};
q = 1;
elseif isnumeric(varargin{end})
x = varargin(2:(end - 1));
n = varargin{end};
V = cell(1,length(x));
q = 1;
elseif ischar(varargin{end})
x = varargin(2:(end - 2));
n = varargin{end - 1};
V = cell(1,length(x));
q = 0;
else
error('improper call')
end
% check transform matrices V = dxdy
%--------------------------------------------------------------------------
for i = 1:length(x)
try
V{i};
catch
V{i} = [];
end
if isempty(V{i}) && any(n == i);
V{i} = speye(spm_length(x{i}));
end
end
% initialise
%--------------------------------------------------------------------------
m = n(end);
xm = spm_vec(x{m});
dx = exp(-8);
J = cell(1,size(V{m},2));
% proceed to derivatives
%==========================================================================
if length(n) == 1
% dfdx
%----------------------------------------------------------------------
f0 = f(x{:});
for i = 1:length(J)
xi = x;
xi{m} = spm_unvec(xm + V{m}(:,i)*dx,x{m});
J{i} = spm_dfdx(f(xi{:}),f0,dx);
end
% return numeric array for first-order derivatives
%======================================================================
% vectorise f
%----------------------------------------------------------------------
f = spm_vec(f0);
% if there are no arguments to differentiate w.r.t. ...
%----------------------------------------------------------------------
if isempty(xm)
J = sparse(length(f),0);
% or there are no arguments to differentiate
%----------------------------------------------------------------------
elseif isempty(f)
J = sparse(0,length(xm));
end
% or differentiation of a vector
%----------------------------------------------------------------------
if isvector(f0) && isnumeric(f0) && q
% concatenate into a matrix
%------------------------------------------------------------------
if size(f0,2) == 1
J = spm_cat(J);
else
J = spm_cat(J')';
end
end
% assign output argument and return
%----------------------------------------------------------------------
varargout{1} = J;
varargout{2} = f0;
else
% dfdxdxdx....
%----------------------------------------------------------------------
f0 = cell(1,length(n));
[f0{:}] = spm_diff(f,x{:},n(1:end - 1),V);
for i = 1:length(J)
xi = x;
xmi = xm + V{m}(:,i)*dx;
xi{m} = spm_unvec(xmi,x{m});
fi = spm_diff(f,xi{:},n(1:end - 1),V);
J{i} = spm_dfdx(fi,f0{1},dx);
end
varargout = [{J} f0];
end
In [ ]:
def spm_dfdx(f,f0,dx):
"""
port of spm_dfdx
...which is defined at the bottom of spm_diff.m
"""
"""
function dfdx = spm_dfdx(f,f0,dx)
% cell subtraction
%--------------------------------------------------------------------------
"""
"""
if iscell(f)
dfdx = f;
for i = 1:length(f(:))
dfdx{i} = spm_dfdx(f{i},f0{i},dx);
end
elseif isstruct(f)
dfdx = (spm_vec(f) - spm_vec(f0))/dx;
else
dfdx = (f - f0)/dx;
end
""";
# Check if input is a dict(cell)
# ...if it is, loop through elements and
# re-call this function
# ...if it isn't,
#
dfdx = []
if type(f) == list:
dfdx = f;
for i = in f[:]:
dfdx{i} = spm_dfdx(f{i},f0{i},dx);
elif isstruct(f)
dfdx = (spm_vec(f) - spm_vec(f0))/dx;
else
dfdx = (f - f0)/dx;
end
""";
In [ ]:
# %load /home/jgriffiths/Code/libraries_of_others/misc/SPM/spm12/spm_diff.m
function [varargout] = spm_diff(varargin)
% matrix high-order numerical differentiation
% FORMAT [dfdx] = spm_diff(f,x,...,n)
% FORMAT [dfdx] = spm_diff(f,x,...,n,V)
% FORMAT [dfdx] = spm_diff(f,x,...,n,'q')
%
% f - [inline] function f(x{1},...)
% x - input argument[s]
% n - arguments to differentiate w.r.t.
%
% V - cell array of matrices that allow for differentiation w.r.t.
% to a linear transformation of the parameters: i.e., returns
%
% df/dy{i}; x = V{i}y{i}; V = dx(i)/dy(i)
%
% q - (char) flag to preclude default concatenation of dfdx
%
% dfdx - df/dx{i} ; n = i
% dfdx{p}...{q} - df/dx{i}dx{j}(q)...dx{k}(p) ; n = [i j ... k]
%
%
% - a cunning recursive routine
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_diff.m 6110 2014-07-21 09:36:13Z karl $
% create inline object
%--------------------------------------------------------------------------
f = spm_funcheck(varargin{1});
% parse input arguments
%--------------------------------------------------------------------------
if iscell(varargin{end})
x = varargin(2:(end - 2));
n = varargin{end - 1};
V = varargin{end};
q = 1;
elseif isnumeric(varargin{end})
x = varargin(2:(end - 1));
n = varargin{end};
V = cell(1,length(x));
q = 1;
elseif ischar(varargin{end})
x = varargin(2:(end - 2));
n = varargin{end - 1};
V = cell(1,length(x));
q = 0;
else
error('improper call')
end
% check transform matrices V = dxdy
%--------------------------------------------------------------------------
for i = 1:length(x)
try
V{i};
catch
V{i} = [];
end
if isempty(V{i}) && any(n == i);
V{i} = speye(spm_length(x{i}));
end
end
% initialise
%--------------------------------------------------------------------------
m = n(end);
xm = spm_vec(x{m});
dx = exp(-8);
J = cell(1,size(V{m},2));
% proceed to derivatives
%==========================================================================
if length(n) == 1
% dfdx
%----------------------------------------------------------------------
f0 = f(x{:});
for i = 1:length(J)
xi = x;
xi{m} = spm_unvec(xm + V{m}(:,i)*dx,x{m});
J{i} = spm_dfdx(f(xi{:}),f0,dx);
end
% return numeric array for first-order derivatives
%======================================================================
% vectorise f
%----------------------------------------------------------------------
f = spm_vec(f0);
% if there are no arguments to differentiate w.r.t. ...
%----------------------------------------------------------------------
if isempty(xm)
J = sparse(length(f),0);
% or there are no arguments to differentiate
%----------------------------------------------------------------------
elseif isempty(f)
J = sparse(0,length(xm));
end
% or differentiation of a vector
%----------------------------------------------------------------------
if isvector(f0) && isnumeric(f0) && q
% concatenate into a matrix
%------------------------------------------------------------------
if size(f0,2) == 1
J = spm_cat(J);
else
J = spm_cat(J')';
end
end
% assign output argument and return
%----------------------------------------------------------------------
varargout{1} = J;
varargout{2} = f0;
else
% dfdxdxdx....
%----------------------------------------------------------------------
f0 = cell(1,length(n));
[f0{:}] = spm_diff(f,x{:},n(1:end - 1),V);
for i = 1:length(J)
xi = x;
xmi = xm + V{m}(:,i)*dx;
xi{m} = spm_unvec(xmi,x{m});
fi = spm_diff(f,xi{:},n(1:end - 1),V);
J{i} = spm_dfdx(fi,f0{1},dx);
end
varargout = [{J} f0];
end
function dfdx = spm_dfdx(f,f0,dx)
% cell subtraction
%--------------------------------------------------------------------------
if iscell(f)
dfdx = f;
for i = 1:length(f(:))
dfdx{i} = spm_dfdx(f{i},f0{i},dx);
end
elseif isstruct(f)
dfdx = (spm_vec(f) - spm_vec(f0))/dx;
else
dfdx = (f - f0)/dx;
end
code:
In [381]:
def spm_vec(X):
# Vectorise a numeric, cell or structure array - a compiled routine
# FORMAT [vX] = spm_vec(X)
#X - numeric, cell or stucture array[s]
#vX - vec(X)
#X = [np.eye(2),3]
# Usage:
#vX = vec(X)
for v_it,v in enumerate(X):
if v_it == 0:
try:
nr,nc = v.shape
vX = np.reshape(v,[nr*nc,1])
except:
vX = np.array([v][:,np.newaxis])
else:
try:
nr,nc = v.shape
vX = np.concatenate([vX,np.reshape(v,[nr*nc,1])])
except:
vX = np.concatenate([vX,np.array([v])[:,np.newaxis]])
return np.squeeze(vX)[:,np.newaxis]
def unvec(Xflat,X):
#nr,nc = X.shape
#vX = np.reshape(X,[nr,nc])
cnt = 0
vXlist = []
for v_it,v in enumerate(X):
try:
nr,nc = v.shape
rng = np.arange(cnt,nr*nc)
vX = np.reshape(Xflat[rng],[nr,nc])
vXlist.append(vX)
cnt+=1
except:
vXlist.append(v)
return vXlist
tests:
matlab spm_vec
In [385]:
%%matlab
X = {eye(2) 3};
vX = spm_vec(X);
disp(['shape: ' num2str(size(vX))]);
vX
python spm_vec
In [386]:
X = [np.eye(2),3]
vX = spm_vec(X)
print 'shape: %s,%s' %vX.shape
print 'value: \n%s' %vX
matlab spm_unvec
In [387]:
%%matlab
X = {eye(2) 3};
vX = spm_unvec(spm_vec(X),X);
vX{:}
python spm_unvec
In [388]:
X = [np.eye(2),3]
vX = spm_unvec(vec(X),X)
#print 'shape: %s,%s' %vX.shape
#print 'value: \n%s' %vX
vX
Out[388]:
code
tests
code
tests
In [ ]:
_,_,D = fx_erp()
In [2]:
%load_ext pymatbridge
In [3]:
%%matlab
addpath(genpath('/home/jgriffiths/Code/libraries_of_others/misc/SPM'));
First half of demo function
In [27]:
%%matlab
rng('default')
Nc = 2; % number of channels
Ns = 2; % number of sources
options.spatial = 'LFP';
options.model = 'ERP';
options.analysis = 'ERP';
M.dipfit.model = options.model;
M.dipfit.type = options.spatial;
M.dipfit.Nc = Nc;
M.dipfit.Ns = Ns;
% sspecify connectivity reciprocal connections with condition specific
% changes in intrinsic and extrinsic connectivity
%--------------------------------------------------------------------------
A{1} = [0 0; 1 0];
A{2} = [0 1; 0 0];
A{3} = [0 0; 0 0];
B{1} = [1 0; 1 0];
C = [1; 0];
[pE,pC] = spm_dcm_neural_priors(A,B,C,options.model);
[gE,gC] = spm_L_priors(M.dipfit);
[x,f] = spm_dcm_x_neural(pE,options.model);
% hyperpriors (assuming a high signal to noise)
%--------------------------------------------------------------------------
hE = 6;
hC = 1/128;
% create model
%--------------------------------------------------------------------------
M.IS = 'spm_gen_erp';
M.G = 'spm_lx_erp';
M.f = f;
M.x = x;
M.pE = pE;
M.pC = pC;
M.gE = gE;
M.gC = gC;
M.hE = hE;
M.hC = hC;
M.m = length(B);
M.n = length(spm_vec(M.x));
M.l = Nc;
M.ns = 64;
% create input structure
%--------------------------------------------------------------------------
dt = 4/1000;
pst = (1:M.ns)*dt;
M.ons = 64;
M.dur = 16;
U.dt = dt;
U.X = [0; 1];
% specified true connectivity (P) and spatial parameters (G) with
% condition specific effects on the intrinsic connectivity of the first
% source and its forward extrinsic connection
%--------------------------------------------------------------------------
P = pE;
G = gE;
P.B{1} = [-1/4 0; 1/2 0];
% generate neuronal response and data
%--------------------------------------------------------------------------
x = spm_gen_erp(P,M,U); % neuronal response
L = spm_lx_erp(G,M.dipfit); % lead field
V = spm_sqrtm(spm_Q(1/2,M.ns)); % square root of noise covariance
for i = 1:length(x)
n = exp(-hE/2)*V*randn(M.ns,Nc); % noise
s{i} = x{i}*L'; % signal
y{i} = s{i} + n; % data (signal plus noise)
end
% data structure specification
%--------------------------------------------------------------------------
Y.y = y;
Y.Q = {spm_Q(1/2,M.ns,1)};
Y.dt = dt;
Y.pst = pst;
Push variables to python
In [31]:
%%matlab -o Y
disp('');
In [34]:
len(Y['y'])
Out[34]:
Run python version
Compare matlab and python values
In [ ]:
In [ ]:
In [ ]:
In [ ]:
a
In [ ]:
In [25]:
%%matlab
% display
%--------------------------------------------------------------------------
spm_figure('Getwin','Figure 1');
subplot(2,1,1)
plot(pst,x{1},'r',pst,x{2},'b')
xlabel('time');ylabel('amplitude');
title('Hidden neuronal states','FontSize',16)
subplot(2,1,2)
plot(pst,s{1},':r',pst,s{2},':b',pst,y{1},'r',pst,y{2},'b')
xlabel('time');ylabel('amplitude');
title('Observed response','FontSize',16)
In [19]:
% Invert model under increasing shrinkage priors on the condition specific
% change in the intrinsic (B) parameter of the first source. This range
% specified by alpha. Because the true value is non-zero, we expect free-
% energy to decrease when the prior covariance falls to 0 and this
% parameter is effectively eliminated. The ensuing estimates are obtained
% while optimising the parameters of the spatial model and the (neuronal)
% extrinsic parameter..
% =========================================================================
In [ ]:
%%matlab
% fix all (neuronal) parameters (except those of interest)
% -------------------------------------------------------------------------
M.pC = spm_unvec(spm_vec(pC)*0,pC);
M.pC.B{1}(1,1) = 1/8;
M.pC.B{1}(2,1) = 1/8;
% vary the prior covariance
% -------------------------------------------------------------------------
alpha = exp(-8:2);
for i = 1:length(alpha)
% reset shrinkage prior
% ---------------------------------------------------------------------
M.pC.B{1}(1,1) = alpha(i);
% full inversion
% ---------------------------------------------------------------------
[Ep,Eg,Cp,Cg,S,F] = spm_nlsi_N(M,U,Y);
Ep_all{i} = spm_vec(Ep);
Cp_all{i} = diag(Cp);
F_all(i) = F;
end
In [ ]:
%%matlab
% Post-hoc reduction for the full (complex) model: this evaluates the free
% energy and posterior distribution over the parameters, given just the
% full prior and posterior for any required reduced prior.
% -------------------------------------------------------------------------
for i = 1:length(alpha)
rC = M.pC;
rC.B{1}(1,1) = alpha(i);
[F,rEp,rCp] = spm_log_evidence_reduce(Ep,Cp,M.pE,M.pC,M.pE,rC);
rEp_all{i} = spm_vec(rEp);
rCp_all{i} = diag(rCp);
rF_all(i) = F;
end
In [12]:
% compare full inversion and model reduction
% -------------------------------------------------------------------------
spm_figure('Getwin','Figure 2');
subplot(2,1,1)
semilogx(alpha,F_all - F_all(end),alpha,rF_all - rF_all(end))
xlabel('prior covariance');ylabel('relative log evidence');
legend('full','posthoc');
title('Full and reduced log evidence','FontSize',16)
% compare parameter estimates
% -------------------------------------------------------------------------
j = find(spm_vec(M.pC));
pP = spm_vec(P);
for i = 1:length(alpha)
fQp(i,:) = Ep_all{i}(j);
fCp(i,:) = Cp_all{i}(j);
rQp(i,:) = rEp_all{i}(j);
vCp(i,:) = rCp_all{i}(j);
Pp(i,:) = pP(j);
end
subplot(2,1,2)
spm_plot_ci(rQp',vCp'), hold on
plot(fQp,'--'), hold on
plot(Pp,'-.'), hold off
xlabel('prior log-covariance');ylabel('difference in log evidence');
title('MAP estimates (full -- reduced - true -.)','FontSize',16)
set(gca,'XTickLabel',log(alpha))
return
% free energy landscape:
% Here, we evaluate the free energy, which is a functional of the data and
% conditional or posterior expectations (noting that posterior precisions
% are functions of the expectations). The free energy can be computed in a
% simple way by inverting the model using a single iteration. In what
% follows, we evaluate the free energy over a range of the two (intrinsic
% and extrinsic) coupling parameters (at the true values of the remaining
% parameters)
% =========================================================================
beta = linspace(-1,1,32); % range of parameters
M.nograph = 1;
M.Nmax = 1;
M.Gmax = 1;
M.Hmax = 1;
M.P = P;
M.Q = G;
for i = 1:length(beta)
for j = 1:length(beta)
% specify posterior expectations (and implicitly precisions)
% -----------------------------------------------------------------
M.P = P;
M.P.B{1}(1,1) = beta(i);
M.P.B{1}(2,1) = beta(j);
% evaluate free energy (without updating expectations)
% -----------------------------------------------------------------
[Ep,Eg,Cp,Cg,S,F] = spm_nlsi_N(M,U,Y);
FF(i,j) = F;
% record conditional covariance at maximum
% -----------------------------------------------------------------
if F >= max(FF(:)), CP = Cp; end
end
end
% apply Occam's window
% -------------------------------------------------------------------------
FF = FF - max(FF(:));
F = max(FF,-64)
% Free energy landscape and associated conditional covariance: the
% conditional covariance is the inverse precision which is proportional
% to the curvature of the variational energy (equivalent to the curvature
% of the free energy functional of posterior or conditional expectations)
% -------------------------------------------------------------------------
spm_figure('Getwin','Figure 3');
j = find(spm_vec(M.pC));
subplot(2,1,1)
imagesc(beta,beta,F), hold on
contour(beta,beta,F,8,'m'), hold on
plot(pP(j(2)),pP(j(1)),'.r','MarkerSize',32), hold off
xlabel('extrinsic');ylabel('intrinsic');
title('Free-energy landscape','FontSize',16)
axis square
subplot(2,1,2)
imagesc(CP(j,j))
xlabel('parameters');ylabel('parameters');
title('Posterior covariance','FontSize',16)
axis square
In [ ]:
%%matlab