Print

Print


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')
>>>         %-#
>>> >
>>>
>>>
>========================================================================