JiscMail Logo
Email discussion lists for the UK Education and Research communities

Help for SPM Archives


SPM Archives

SPM Archives


SPM@JISCMAIL.AC.UK


View:

Message:

[

First

|

Previous

|

Next

|

Last

]

By Topic:

[

First

|

Previous

|

Next

|

Last

]

By Author:

[

First

|

Previous

|

Next

|

Last

]

Font:

Proportional Font

LISTSERV Archives

LISTSERV Archives

SPM Home

SPM Home

SPM  February 2007

SPM February 2007

Options

Subscribe or Unsubscribe

Subscribe or Unsubscribe

Log In

Log In

Get Password

Get Password

Subject:

Re: spm5 bayesian problems

From:

Martin Kronbichler <[log in to unmask]>

Reply-To:

Martin Kronbichler <[log in to unmask]>

Date:

Wed, 7 Feb 2007 09:47:22 +0000

Content-Type:

text/plain

Parts/Attachments:

Parts/Attachments

text/plain (701 lines)

I have playing around with my data set a little more, interestingly
the error (detailed below) does occur when is specifiy the contrast
after the bayesian estimation but not when i specify the contrast before
estimation (i.e. when the contrast is calculated during the estimation
stage) altough i understand that the general recommendation for bayesian
fmri-analysis is to specify contrasts before estimating the model i'm still
puzzled that calculating contrasts after estimation can lead to such odd and
invalid results? Is this a bug or are such errors expected when one
calculates the contrast after the bayesian estimation (my problem is that
i'm not always sure which contrasts i want to test and it the bayesian
single-subject estimation takes quite some time). 
Another question: since i do not seem to be able to specify an effect size
threshold for contrast which are calculated during the estimation stage, how
is the effect size cut-off determind for such contrast

greetings,

martin


On Wed, 31 Jan 2007 13:48:44 +0000, Martin Kronbichler
<[log in to unmask]> wrote:

>Hi Darren,
>
>no the PPMs (all brain voxels) are still black whatever value i choose for 
>the effectsize threshold and/or the posterior proability. 
>
>On Wed, 31 Jan 2007 07:33:50 -0600, d gitelman <d-
>[log in to unmask]> wrote:
>
>>Dear Martin
>>
>>This is just a warning, not an error. Do the maps look correct now- i.e.,
>>not totally black and the thresholds appear reasonable. If so then it's
>>fine.
>>
>>darren
>>
>>> -----Original Message-----
>>> From: Martin Kronbichler [mailto:[log in to unmask]]
>>> Sent: Wednesday, January 31, 2007 4:18 AM
>>> To: Darren G
>>> Cc: Martin Kronbichler
>>> Subject: Re: spm5 bayesian problems
>>>
>>> Dear Darren,
>>>
>>> thanks for your fast response. However, the fix didn't work for me.
>>> I tried to reestimate the design with you updated
>>> spm_spm_Bayes as well as manually replacing "tmp = tmp +
>>> NaN*(~tmp); " with "tmp(~tmp) = NaN;" in both lines in my
>>> local spm_spm_Bayes. In both cases i get the same error when
>>> i test a specific contrast
>>> (3 conditions -1 0.5 0.5), however when i use a different
>>> contrast (e.g., condition 1 vs. implicit baseline: 1 0 0) i
>>> don't get the error mesage and the results look reasonable.
>>> It seems as if the error occurs whenever i try to contrast
>>> any of the modelled conditions against one of the other
>>> modelled conditions, but not when i estimate the design,
>>> Strangely, the same design gives reasonable results in an
>>> other patients dataset and normally first-level bayesian
>>> analysis works extremly well for my data-sets. You can see
>>> the error returnded when i use dbstop below, to my limited
>>> understanding the error starts in spm_Ncdf (the cum normal
>>> distribution..). I might add that i use (still) MATLAB R13).
>>> Maybe you have any ideas left what's going wrong?
>>>
>>>
>>> Greetings,
>>>
>>> Martin
>>>
>>>
>>>
>>> In /opt/spm5/spm_Ncdf.m at line 80
>>>   In /opt/spm5/spm_contrasts.m at line 214
>>>   In /opt/spm5/spm_getSPM.m at line 474
>>>   In /opt/spm5/spm_results_ui.m at line 264
>>>   In /opt/spm5/spm_input.m at line 1121
>>>   In /opt/spm5/spm_getSPM.m at line 376
>>>   In /opt/spm5/spm_results_ui.m at line 264
>>>
>>> Warning in ==> /opt/spm5/spm_Ncdf.m
>>> On line 80  ==>         warning('Returning NaN for out of
>>> range arguments'), end
>>>
>>>
>>> On Tue, 30 Jan 2007 10:03:24 -0600, d gitelman
>>> <[log in to unmask]>
>>> wrote:
>>>
>>> > Dear Martin
>>> >
>>> >Greetings. The problem is in spm_spm_Bayes.m when one uses
>>> Matlab R14
>>> >or greater. If you go to approx lines 294 and 301 this code :
>>> >     tmp       = tmp + NaN*(~tmp);
>>> >
>>> >produces the NaN's.
>>> >I have updated this in the FIL version of the file and hopefully it
>>> >will make its way to wider distribution at a future update
>>> (or it may
>>> >have already).
>>> >
>>> >In any case replace the above with this line in both places.
>>> >
>>> >    tmp(~tmp) = NaN;
>>> >
>>> >You will have to re-estimate the dataset, but then all
>>> should be fine.
>>> >
>>> >Attached is also an updated version of spm_spm_Bayes.m
>>> >
>>> >Regards,
>>> >Darren
>>> >
>>> >> -----Original Message-----
>>> >> From: SPM (Statistical Parametric Mapping)
>>> >> [mailto:[log in to unmask]] On Behalf Of Martin Kronbichler
>>> >> Sent: Tuesday, January 30, 2007 4:12 AM
>>> >> To: [log in to unmask]
>>> >> Subject: [SPM] spm5 bayesian problems
>>> >>
>>> >> Dear SPM Experts,
>>> >>
>>> >> i have recently started to use variational bayesian first level
>>> >> analyses for analysing presurgical language and motor experiments
>>> >> because i had the feeling that the posterior probalities
>>> are easier
>>> >> to communicate and easier to understand (especially when making
>>> >> inferences on how likely increased activatity for speech
>>> or motor is
>>> >> in regions near to the lesions).
>>> >> Generally, the bayesian analyses worked extremly well.
>>> >> However i have stumpled across a strange result i do not
>>> understand
>>> >> when analysing a block-design language experiment
>>> >> (2 blockedc onditions: simple tones vs. auditory words
>>> presented in 2
>>> >> sessions). Whatever effect size threshold (even completerly
>>> >> unrealistic values like > 1000) or posterior probality i
>>> choose the
>>> >> PPM shows the whole brains as exhibiting posterior prob.
>>> of 1 (i.e.
>>> >> the whole brain is "activated"). SPM5 displays the following error
>>> >> message when is estimate the contrasts:
>>> >>
>>> >> spm{P} image  1                 :                   ...computing
>>> >> Warning: Returning NaN for out of range arguments
>>> >> > In /opt/spm5/spm_Ncdf.m at line 80
>>> >>   In /opt/spm5/spm_contrasts.m at line 214
>>> >>   In /opt/spm5/spm_getSPM.m at line 474
>>> >>   In /opt/spm5/spm_results_ui.m at line 264
>>> >>       ...written spmP_0001.img
>>> >>         SPM computation                 :
>>> >> ...initialising
>>> >>         SPM computation                 :
>>> >>    ...done
>>> >>         SPM computation                 :
>>> >>    ...done
>>> >>
>>> >> I have tried to remove and add motion parameters, analyse only one
>>> >> session, use different basis functions etc... but whatever
>>> i do the
>>> >> problem persits.
>>> >> Inspection of the files generated by SPM5 shows that the beta and
>>> >> SDbeta images look reasonable. The spmP and con_SD images however
>>> >> contain only NaNs..
>>> >> I might add that a conventional analysis results in stastically
>>> >> highly signficant voxels...
>>> >> Anybody has any idea what is causing this problem?
>>> >>
>>> >> Greetings from Salzburg,
>>> >>
>>> >> Martin
>>> >>
>>> >>
>>> >> Martin Kronbichler,  M.Sc.
>>> >> ---------------------------------------------
>>> >> Department of Psychology
>>> >>
>>> >>
>>> >> Center for Neurocognitive Research
>>> >> University of Salzburg
>>> >>
>>> >>     Hellbrunnerstr.34, A 5020 Salzburg, Austria
>>> >> e-mail:[log in to unmask]
>>> >> Tel.:+43/(0)/662/8044-5162
>>> >>
>>> >> Fax:+43/(0)/662/8044-5126
>>> >>
>>> >> ---------------------------------------------
>>> >> Department of Neurology
>>> >> Center for Neurocognitive Research  Christian-Doppler Clinic,
>>> >> Paracelsus Private Medical University Ignaz Harrerstr.79, A 5020
>>> >> Salzburg, Austria e-mail:[log in to unmask]
>>> >> Tel.:+43/(0)/662/4483-3966
>>> >> ---------------------------------------------
>>> >>
>>> >
>>> >function [SPM] = spm_spm_Bayes(SPM)
>>> >% Conditional parameter estimation of a General Linear Model
>>> % FORMAT
>>> >[SPM] = spm_spm_Bayes(SPM)
>>> >%____________________________________________________________
>>> __________
>>> >_
>>> >%
>>> >% For single subject fMRI analysis this function switches to % using
>>> >voxel-wise GLM-AR models that are spatially regularised %
>>> using the VB
>>> >framework. This is implemented using spm_spm_vb.m.
>>> >% One need not have previously analysed the data using
>>> classical estimation.
>>> >% For more on this option please see the help in that file.
>>> >% Otherwise, read on.
>>> >%
>>> >% spm_spm_Bayes returns to voxels identified by spm_spm (ML
>>> parameter %
>>> >estimation) to get conditional parameter estimates and ReML hyper- %
>>> >parameter estimates.  These estimates use prior covariances,
>>> on the %
>>> >parameters, from emprical Bayes.  These PEB prior variances
>>> come from %
>>> >the hierarchical model that obtains by considering voxels as
>>> providing
>>> >% a second level.  Put simply, the variance in parameters,
>>> over voxels,
>>> >% is used as a prior variance from the point of view of any
>>> one voxel.
>>> >% The error covariance hyperparameters are re-estimated in
>>> the light of
>>> >% these priors.  The approach adopted is essentially a fully
>>> Bayesian %
>>> >analysis at each voxel, using emprical Bayesian prior variance %
>>> >estimators over voxels.
>>> >%
>>> >% Each separable partition (i.e. session) is assigned its own %
>>> >hyperparameter but within session covariance components are lumped %
>>> >together, using their relative expectations over voxels.
>>> This makes %
>>> >things much more computationally efficient and avoids inefficient %
>>> >voxel-specific multiple hyperparameter estimates.
>>> >%
>>> >% spm_spm_Bayes adds the following fields to SPM:
>>> >%
>>> >%                           ----------------
>>> >%
>>> >%
>>> >%	SPM.PPM.l      = session-specific hyperparameter means
>>> >%	SPM.PPM.Cb     = empirical prior parameter covariances
>>> >%	SPM.PPM.C      = conditional covariances of parameters
>>> >%	SPM.PPM.dC{i}  = dC/dl;
>>> >%	SPM.PPM.ddC{i} = ddC/dldl
>>> >%
>>> >% The derivatives are used to compute the conditional variance of
>>> >various % contrasts in spm_getSPM, using a Taylor expansion
>>> about the
>>> >hyperparameter % means.
>>> >%
>>> >%
>>> >%                           ----------------
>>> >%
>>> >%	SPM.VCbeta     - Handles of conditional parameter estimates
>>> >%	SPM.VHp        - Handles of hyperparameter estimates
>>> >%
>>> >%                           ----------------
>>> >%
>>> >% Cbeta_????.{img,hdr}                     - conditional
>>> parameter images
>>> >% These are 16-bit (float) images of the conditional estimates. The
>>> >image % files are numbered according to the corresponding
>>> column of the
>>> >% design matrix. Voxels outside the analysis mask (mask.img)
>>> are given
>>> >% value NaN.
>>> >%
>>> >%                           ----------------
>>> >%
>>> >% CHp_????.{img,hdr}              - error covariance
>>> hyperparamter images
>>> >% This is a 32-bit (double) image of the ReML error variance
>>> estimate.
>>> >% for each separable partition (Session).  Voxels outside
>>> the analysis
>>> >% mask are given value NaN.
>>> >%
>>> >%____________________________________________________________
>>> __________
>>> >_ % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
>>> >
>>> >% Karl Friston
>>> >% $Id: spm_spm_Bayes.m 587 2006-08-07 04:38:22Z Darren $
>>> >
>>> >
>>> >%-Say hello
>>> >%------------------------------------------------------------
>>> ----------
>>> >- Finter = spm('FigName','Stats: Bayesian estimation...');
>>> >
>>> >%-Select SPM.mat & change directory
>>> >%------------------------------------------------------------
>>> ----------
>>> >-
>>> >if ~nargin
>>> >	swd = spm_str_manip(spm_select(1,'^SPM\.mat$','Select
>>> SPM.mat'),'H');
>>> >	load(fullfile(swd,'SPM.mat'))
>>> >	cd(swd)
>>> >end
>>> >
>>> >try
>>> >	M      = SPM.xVol.M;
>>> >	DIM    = SPM.xVol.DIM;
>>> >	xdim   = DIM(1); ydim = DIM(2); zdim = DIM(3);
>>> >	XYZ    = SPM.xVol.XYZ;
>>> >catch
>>> >	helpdlg({	'Please do a ML estimation first',...
>>> >			'This identifies the voxels to analyse'});
>>> >	spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
>>> >	return
>>> >end
>>> >
>>> >
>>> >%============================================================
>>> ===========
>>> >% - A N A L Y S I S   P R E L I M I N A R I E S
>>> >%============================================================
>>> ==========
>>> >=
>>> >
>>> >%-Initialise output images
>>> >%============================================================
>>> ===========
>>> >fprintf('%-40s: %30s','Output images','...initialising')
>>>         %-#
>>> >
>>> >%-Intialise oonditional estimate image files
>>> >%------------------------------------------------------------
>>> -----------
>>> >xX             = SPM.xX;
>>> >[nScan nBeta]  = size(xX.X);
>>> >Vbeta(1:nBeta) = deal(struct(...
>>> >			'fname',	[],...
>>> >			'dim',		DIM',...
>>> >			'dt',		[spm_type('float32'),
>>> spm_platform('bigend')],...
>>> >			'mat',		M,...
>>> >			'pinfo',	[1 0 0]',...
>>> >			'descrip',	''));
>>> >for i = 1:nBeta
>>> >	Vbeta(i).fname   = sprintf('Cbeta_%04d.img',i);
>>> >	Vbeta(i).descrip = sprintf('Cond. beta (%04d) -
>>> %s',i,xX.name{i});
>>> >	spm_unlink(Vbeta(i).fname)
>>> >end
>>> >Vbeta = spm_create_vol(Vbeta);
>>> >
>>> >%-Intialise ReML hyperparameter image files
>>> >%------------------------------------------------------------
>>> ----------
>>> >-
>>> >try
>>> >	nHp       = length(SPM.nscan);
>>> >catch
>>> >	nHp       = nScan;
>>> >	SPM.nscan = nScan;
>>> >end
>>> >
>>> >VHp(1:nHp)        = deal(struct(...
>>> >			'fname',	[],...
>>> >			'dim',		DIM',...
>>> >			'dt',		[spm_type('float64'),
>>> spm_platform('bigend')],...
>>> >			'mat',		M,...
>>> >			'pinfo',	[1 0 0]',...
>>> >			'descrip',	''));
>>> >for i = 1:nHp
>>> >	VHp(i).fname   = sprintf('Hp_%04d.img',i);
>>> >	VHp(i).descrip = sprintf('Hyperparameter (%04d)',i);
>>> >	spm_unlink(VHp(i).fname)
>>> >end
>>> >VHp   = spm_create_vol(VHp);
>>> >
>>> >fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...initialised
>>')        %-#
>>> >
>>> >
>>> >%============================================================
>>> ==========
>>> >= % - E M P I R I C A L  B A Y E S  F O R  P R I O R  V A R
>>> I A N C E
>>> >%============================================================
>>> ===========
>>> >fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...estimating
>>> priors')   %-#
>>> >
>>> >% get row u{i} and column v{i}/v0{i} indices for separable designs
>>> >%------------------------------------------------------------
>>> ----------
>>> >s      = nHp;
>>> >if isfield(SPM,'Sess')
>>> >	for i = 1:s
>>> >		 u{i} = SPM.Sess(i).row;
>>> >		 v{i} = SPM.Sess(i).col;
>>> >		v0{i} = xX.iB(i);
>>> >	end
>>> >else
>>> >	 u{1} = [1:nScan];
>>> >	 v{1} = [xX.iH xX.iC];
>>> >	v0{1} = [xX.iB xX.iG];
>>> >end
>>> >
>>> >% cycle over separarable partitions
>>> >%------------------------------------------------------------
>>> ----------
>>> >-
>>> >for i = 1:s
>>> >
>>> >	% Get design X and confounds X0
>>> >	%---------------------------------------------------------------
>>> >	fprintf('%-30s- %i\n','  ReML Session',i);
>>> >	X     = xX.X(u{i}, v{i});
>>> >	X0    = xX.X(u{i},v0{i});
>>> >	[m n] = size(X);
>>> >
>>> >	% add confound in 'filter'
>>> >	%---------------------------------------------------------------
>>> >	if isstruct(xX.K)
>>> >		X0 = full([X0 xX.K(i).X0]);
>>> >	end
>>> >
>>> >	% orthogonalize X w.r.t. X0
>>> >	%---------------------------------------------------------------
>>> >	X     = X - X0*(pinv(X0)*X);
>>> >
>>> >	% covariance components induced by parameter variations {Q}
>>> >	%---------------------------------------------------------------
>>> >	for j = 1:n
>>> >		Q{j} = X*sparse(j,j,1,n,n)*X';
>>> >	end
>>> >
>>> >	% covariance components induced by error non-sphericity {V}
>>> >	%---------------------------------------------------------------
>>> >	Q{n + 1} = SPM.xVi.V(u{i},u{i});
>>> >
>>> >	% ReML covariance component estimation
>>> >	%---------------------------------------------------------------
>>> >	[C h]   = spm_reml(SPM.xVi.CY,X0,Q);
>>> >
>>> >	% check for negative variance components
>>> >	%---------------------------------------------------------------
>>> >	h       = abs(h);
>>> >
>>> >	% 2-level model for this partition using prior variances sP(i)
>>> >	% treat confounds as fixed (i.e. infinite prior variance)
>>> >	%---------------------------------------------------------------
>>> >	n0      = size(X0,2);
>>> >	Cb      = blkdiag(diag(h(1:n)),speye(n0,n0)*1e8);
>>> >	P{1}.X  = [X X0];
>>> >	P{1}.C  = {SPM.xVi.V};
>>> >	P{2}.X  = sparse(size(P{1}.X,2),1);
>>> >	P{2}.C  = Cb;
>>> >
>>> >	sP(i).P = P;
>>> >	sP(i).u = u{:};
>>> >	sP(i).v = [v{:} v0{:}];
>>> >end
>>> >
>>> >
>>> >%============================================================
>>> ===========
>>> >% - F I T   M O D E L   &   W R I T E   P A R A M E T E R
>>> I M A G E S
>>> >%============================================================
>>> ==========
>>> >=
>>> >
>>> >%-Cycle to avoid memory problems (plane by plane)
>>> >%============================================================
>>> ==========
>>> >= spm_progress_bar('Init',100,'Bayesian estimation','');
>>> helpdlg('This
>>> >could take some time');
>>> >spm('Pointer','Watch')
>>> >
>>> >%-maxMem is the maximum amount of data processed at a time (bytes)
>>> >%------------------------------------------------------------
>>> ----------
>>> >-
>>> >global defaults
>>> >MAXMEM = defaults.stats.maxmem;
>>> >blksz  = ceil(MAXMEM/8/nScan);
>>> >SHp    = 0;				% sum of hyperparameters
>>> >for  z = 1:zdim
>>> >
>>> >    % current plane-specific parameters
>>> >
>>> %-------------------------------------------------------------------
>>> >    U       = find(XYZ(3,:) == z);
>>> >    nbch    = ceil(length(U)/blksz);
>>> >    CrBl    = zeros(nBeta,length(U));	%-conditional
>>> parameter estimates
>>> >    CrHp    = zeros(nHp,  length(U));	%-ReML
>>> hyperparameter estimates
>>> >    for bch = 1:nbch			%-loop over bunches of
>>> lines (planks)
>>> >
>>> >	%-construct list of voxels in this block
>>> >	%---------------------------------------------------------------
>>> >	I     = [1:blksz] + (bch - 1)*blksz;
>>> >	I     = I(I <= length(U));
>>> >	xyz   = XYZ(:,U(I));
>>> >	nVox  = size(xyz,2);
>>> >
>>> >	%-Get response variable
>>> >	%---------------------------------------------------------------
>>> >	Y     = spm_get_data(SPM.xY.VY,xyz);
>>> >
>>> >	%-Conditional estimates (per partition, per voxel)
>>> >	%---------------------------------------------------------------
>>> >	beta  = zeros(nBeta,nVox);
>>> >	Hp    = zeros(nHp,  nVox);
>>> >	for j = 1:s
>>> >		P     = sP(j).P;
>>> >		u     = sP(j).u;
>>> >		v     = sP(j).v;
>>> >		for i = 1:nVox
>>> >			[C P]      = spm_PEB(Y(u,i),P);
>>> >			beta(v,i)  = C{2}.E(1:length(v));
>>> >			Hp(j,i)    = C{1}.h;
>>> >		end
>>> >	end
>>> >
>>> >	%-Save for current plane in memory as we go along
>>> >	%---------------------------------------------------------------
>>> >	CrBl(:,I) = beta;
>>> >	CrHp(:,I) = Hp;
>>> >	SHp       = SHp + sum(Hp,2);
>>> >
>>> >    end % (bch)
>>> >
>>> >
>>> >    %-write out plane data to image files
>>> >
>>> > %===================================================================
>>> >
>>> >    %-Write conditional beta images
>>> >
>>> %-------------------------------------------------------------------
>>> >    for i = 1:nBeta
>>> >	tmp       = sparse(XYZ(1,U),XYZ(2,U),CrBl(i,:),xdim,ydim);
>>> >%     The following line (and the same code on line 306
>>> below) results in
>>> >%     images with all NaN's (in Matlab R14  or greater). It has been
>>> >%     commented out by DRG, and replaced by the line:
>>> tmp(~tmp) = NaN;
>>> >%     tmp       = tmp + NaN*(~tmp);
>>> >    tmp(~tmp) = NaN;
>>> >	Vbeta(i)  = spm_write_plane(Vbeta(i),tmp,z);
>>> >    end
>>> >
>>> >    %-Write hyperparameter images
>>> >
>>> %-------------------------------------------------------------------
>>> >    for i = 1:nHp
>>> >	tmp       = sparse(XYZ(1,U),XYZ(2,U),CrHp(i,:),xdim,ydim);
>>> >%         tmp       = tmp + NaN*(~tmp);
>>> >        tmp(~tmp) = NaN;
>>> >	VHp(i)    = spm_write_plane(VHp(i),tmp,z);
>>> >    end
>>> >
>>> >
>>> >    %-Report progress
>>> >
>>> %-------------------------------------------------------------------
>>> >    spm_progress_bar('Set',100*(z - 1)/zdim);
>>> >
>>> >
>>> >end % (for z = 1:zdim)
>>> >fprintf('\n')
>>>         %-#
>>> >spm_progress_bar('Clear')
>>> >
>>> >%============================================================
>>> ===========
>>> >% - P O S T   E S T I M A T I O N
>>> >%============================================================
>>> ==========
>>> >=
>>> >
>>> >% Taylor expansion for conditional covariance
>>> >%------------------------------------------------------------
>>> ----------
>>> >-
>>> >fprintf('%-40s: %30s\n','Non-sphericity','...REML estimation') %-#
>>> >
>>> >% expansion point (mean hyperparameters)
>>> >%------------------------------------------------------------
>>> -----------
>>> >l     = SHp/SPM.xVol.S;
>>> >
>>> >% change in conditional coavriance w.r.t. hyperparameters
>>> >%------------------------------------------------------------
>>> -----------
>>> >n     = size(xX.X,2);
>>> >PPM.l = l;
>>> >for i = 1:s
>>> >	PPM.dC{i}  = sparse(n,n);
>>> >	PPM.ddC{i} = sparse(n,n);
>>> >end
>>> >for i = 1:s
>>> >
>>> >	P     = sP(i).P;
>>> >	u     = sP(i).u;
>>> >	v     = sP(i).v;
>>> >
>>> >	% derivatives of conditional covariance w.r.t. hyperparameters
>>> >	%---------------------------------------------------------------
>>> >    if spm_matlab_version_chk('6.5.1') < 0
>>> >        % MATLAB 6.5.0 R13 has a bug in dealing with sparse
>>> matrices so
>>> revert to full matrices
>>> >        fprintf('%-40s','MATLAB 6.5.0 R13: must use full matrices -
>>> > this
>>> could cause memory problems...');
>>> >        d     = full(P{1}.X'*inv(P{1}.C{1})*P{1}.X);
>>> >        Cby   = full(inv(d/l(i) + inv(P{2}.C)));
>>> >        d     = full(d*Cby);
>>> >        dC    = full(Cby*d/(l(i)^2));
>>> >        ddC   = full(2*(dC/(l(i)^2) - Cby/(l(i)^3))*d);
>>> >    else
>>> >        % all other MATLAB versions: use sparse matrices
>>> >        d     = P{1}.X'*inv(P{1}.C{1})*P{1}.X;
>>> >        Cby   = inv(d/l(i) + inv(P{2}.C));
>>> >        d     = d*Cby;
>>> >        dC    = Cby*d/(l(i)^2);
>>> >        ddC   = 2*(dC/(l(i)^2) - Cby/(l(i)^3))*d;
>>> >    end
>>> >
>>> >	% place in output structure
>>> >	%---------------------------------------------------------------
>>> >    if spm_matlab_version_chk('6.5.1') < 0,
>>> >        % MATLAB 6.5.0 R13: revert to full matrices
>>> >        j               = 1:length(v);
>>> >        PPM.Cb(v,v)     = full(P{2}.C(j,j));
>>> >        PPM.Cby(v,v)    = full(Cby(j,j));
>>> >        PPM.dC{i}(v,v)  = full(dC(j,j));
>>> >        PPM.ddC{i}(v,v) = full(ddC(j,j));
>>> >    else
>>> >        % all other MATLAB versions: use sparse matrices
>>> >        j               = 1:length(v);
>>> >        PPM.Cb(v,v)     = P{2}.C(j,j);
>>> >        PPM.Cby(v,v)    = Cby(j,j);
>>> >        PPM.dC{i}(v,v)  = dC(j,j);
>>> >        PPM.ddC{i}(v,v) = ddC(j,j);
>>> >    end
>>> >
>>> >end
>>> >
>>> >%-Save remaining results files and analysis parameters
>>> >%============================================================
>>> ===========
>>> >fprintf('%-40s: %30s','Saving results','...writing')
>>>         %-#
>>> >
>>> >%-Save analysis parameters in SPM.mat file
>>> >%------------------------------------------------------------
>>> -----------
>>> >SPM.VCbeta = Vbeta;			% Filenames - parameters
>>> >SPM.VHp    = VHp;			% Filenames - hyperparameters
>>> >SPM.PPM    = PPM;			% PPM structure
>>> >
>>> >if spm_matlab_version_chk('7') >=0
>>> >	save('SPM', 'SPM', '-V6');
>>> >else
>>> >	save('SPM', 'SPM');
>>> >end;
>>> >
>>> >fprintf('%s%30s\n',repmat(sprintf('\b'),1,30),'...done')
>>>           %-#
>>> >
>>> >
>>> >%============================================================
>>> ==========
>>> >=
>>> >%- E N D: Cleanup GUI
>>> >%============================================================
>>> ==========
>>> >=
>>> >spm('FigName','Stats: done',Finter); spm('Pointer','Arrow')
>>> >fprintf('%-40s: %30s\n','Completed',spm('time'))
>>>         %-#
>>> >fprintf('...use the results section for assessment\n\n')
>>>         %-#
>>> >
>>>
>>>
>========================================================================

Top of Message | Previous Page | Permalink

JiscMail Tools


RSS Feeds and Sharing


Advanced Options


Archives

March 2024
February 2024
January 2024
December 2023
November 2023
October 2023
September 2023
August 2023
July 2023
June 2023
May 2023
April 2023
March 2023
February 2023
January 2023
December 2022
November 2022
October 2022
September 2022
August 2022
July 2022
June 2022
May 2022
April 2022
March 2022
February 2022
January 2022
December 2021
November 2021
October 2021
September 2021
August 2021
July 2021
June 2021
May 2021
April 2021
March 2021
February 2021
January 2021
December 2020
November 2020
October 2020
September 2020
August 2020
July 2020
June 2020
May 2020
April 2020
March 2020
February 2020
January 2020
December 2019
November 2019
October 2019
September 2019
August 2019
July 2019
June 2019
May 2019
April 2019
March 2019
February 2019
January 2019
December 2018
November 2018
October 2018
September 2018
August 2018
July 2018
June 2018
May 2018
April 2018
March 2018
February 2018
January 2018
December 2017
November 2017
October 2017
September 2017
August 2017
July 2017
June 2017
May 2017
April 2017
March 2017
February 2017
January 2017
December 2016
November 2016
October 2016
September 2016
August 2016
July 2016
June 2016
May 2016
April 2016
March 2016
February 2016
January 2016
December 2015
November 2015
October 2015
September 2015
August 2015
July 2015
June 2015
May 2015
April 2015
March 2015
February 2015
January 2015
December 2014
November 2014
October 2014
September 2014
August 2014
July 2014
June 2014
May 2014
April 2014
March 2014
February 2014
January 2014
December 2013
November 2013
October 2013
September 2013
August 2013
July 2013
June 2013
May 2013
April 2013
March 2013
February 2013
January 2013
December 2012
November 2012
October 2012
September 2012
August 2012
July 2012
June 2012
May 2012
April 2012
March 2012
February 2012
January 2012
December 2011
November 2011
October 2011
September 2011
August 2011
July 2011
June 2011
May 2011
April 2011
March 2011
February 2011
January 2011
December 2010
November 2010
October 2010
September 2010
August 2010
July 2010
June 2010
May 2010
April 2010
March 2010
February 2010
January 2010
December 2009
November 2009
October 2009
September 2009
August 2009
July 2009
June 2009
May 2009
April 2009
March 2009
February 2009
January 2009
December 2008
November 2008
October 2008
September 2008
August 2008
July 2008
June 2008
May 2008
April 2008
March 2008
February 2008
January 2008
December 2007
November 2007
October 2007
September 2007
August 2007
July 2007
June 2007
May 2007
April 2007
March 2007
February 2007
January 2007
2006
2005
2004
2003
2002
2001
2000
1999
1998


JiscMail is a Jisc service.

View our service policies at https://www.jiscmail.ac.uk/policyandsecurity/ and Jisc's privacy policy at https://www.jisc.ac.uk/website/privacy-notice

For help and support help@jisc.ac.uk

Secured by F-Secure Anti-Virus CataList Email List Search Powered by the LISTSERV Email List Manager