pr_fmristat_ar

 function returns estimated AR coefficients using fmristat algorithm
 FORMAT [rho,Vmhalf,V] = pr_fmristat_ar(res,X,nlags)
 
 See http://www.math.mcgill.ca/keith/fmristat/ and
 fmrilm.m in fmristat package for code, and
 Worsley, K.J., Liao, C., Aston, J., Petre, V., Duncan, G.H., Morales,
 F., Evans, A.C. (2002). A general statistical analysis for fMRI
 data. NeuroImage, 15:1-15 - for description of the algorithm
  
 $Id$

0001 function [rho,Vmhalf,V] = pr_fmristat_ar(res,X,nlags)
0002 % function returns estimated AR coefficients using fmristat algorithm
0003 % FORMAT [rho,Vmhalf,V] = pr_fmristat_ar(res,X,nlags)
0004 %
0005 % See http://www.math.mcgill.ca/keith/fmristat/ and
0006 % fmrilm.m in fmristat package for code, and
0007 % Worsley, K.J., Liao, C., Aston, J., Petre, V., Duncan, G.H., Morales,
0008 % F., Evans, A.C. (2002). A general statistical analysis for fMRI
0009 % data. NeuroImage, 15:1-15 - for description of the algorithm
0010 %
0011 % $Id$
0012 
0013 % This is the copyright notice from fmrilm:
0014 %############################################################################
0015 % COPYRIGHT:   Copyright 2002 K.J. Worsley
0016 %              Department of Mathematics and Statistics,
0017 %              McConnell Brain Imaging Center,
0018 %              Montreal Neurological Institute,
0019 %              McGill University, Montreal, Quebec, Canada.
0020 %              worsley@math.mcgill.ca, liao@math.mcgill.ca
0021 %
0022 %              Permission to use, copy, modify, and distribute this
0023 %              software and its documentation for any purpose and without
0024 %              fee is hereby granted, provided that the above copyright
0025 %              notice appear in all copies.  The author and McGill University
0026 %              make no representations about the suitability of this
0027 %              software for any purpose.  It is provided "as is" without
0028 %              express or implied warranty.
0029 %############################################################################
0030   
0031 if nargin < 2
0032   error('Need covariance and design');
0033 end
0034 if nargin < 3
0035   nlags = 1;
0036 end
0037 
0038 sX = spm_sp('Set', X);
0039 R  = spm_sp('r', sX);
0040 
0041 nlp1     = nlags+1;
0042 [n nvox] = size(res);
0043 
0044 % Bias reduction
0045 M=zeros(nlp1);
0046 for i=1:(nlp1)
0047   for j=1:(nlp1)
0048     Di=(diag(ones(1,n-i+1),i-1)+diag(ones(1,n-i+1),-i+1))/(1+(i==1));
0049     Dj=(diag(ones(1,n-j+1),j-1)+diag(ones(1,n-j+1),-j+1))/(1+(j==1));
0050     M(i,j)=trace(R*Di*R*Dj)/(1+(i>1));
0051   end
0052 end
0053 invM = inv(M);
0054 
0055 a = zeros(nlp1,nvox);
0056 for lag = 0:nlags
0057   a(lag+1,:)= sum(res(1:(n-lag),:).*res((lag+1):n,:));
0058 end
0059 
0060 vhat = invM*a;
0061 rho  = vhat(2:end,:) ./ (ones(nlags, 1) * vhat(1,:));
0062 rho  = mean(rho,2)';
0063 
0064 if nargout > 1
0065   % Whitening matrix; Appendix A3 Worsley et al (2002)
0066   % Modified according to fmrilm code
0067   [Ainvt posdef] = chol(toeplitz([1 rho]));
0068   nl=size(Ainvt,1);
0069   A=inv(Ainvt');
0070   Vmhalf = zeros(n,n);
0071   B=ones(n-nl,1)*A(nl,:);
0072   Vmhalf(nl+1:end,:) = spdiags(B,1:nl,n-nl,n);
0073   Vmhalf(1:nl,1:nl)  = A;
0074 end
0075 
0076 if nargout > 2
0077   % Estimated covariance
0078   V = inv(Vmhalf*Vmhalf');
0079 end