Hi Martin:
Some answers:
the threshold is specified at the contrast definition stage. The default is
1 conditional SD of the prior covariance, though you can choose other
thresholds.
I still don't understand why you are having trouble before and after
estimation in specifying contrasts, or why the images should still be black.
I will email you my updated files dealing with Bayes off list to let you
check if they work. If so I can then update the list.
Darren
> -----Original Message-----
> From: SPM (Statistical Parametric Mapping)
> [mailto:[log in to unmask]] On Behalf Of Martin Kronbichler
> Sent: Wednesday, February 07, 2007 3:47 AM
> To: [log in to unmask]
> Subject: Re: [SPM] spm5 bayesian problems
>
> 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')
> >>> %-#
> >>> >
> >>>
> >>>
> >=============================================================
> ==========
> >=
>
