spm_vb_init_volume

 Initialise generic aspects of slice structure for VB GLMAR models
 FORMAT [slice] = spm_vb_init_volume (X,p)

 X design matrix
 p AR model order

 %W% Will Penny and Nelson Trujillo-Barreto %E%

0001  function [slice] = spm_vb_init_volume (X,p)
0002 % Initialise generic aspects of slice structure for VB GLMAR models
0003 % FORMAT [slice] = spm_vb_init_volume (X,p)
0004 %
0005 % X design matrix
0006 % p AR model order
0007 %
0008 % %W% Will Penny and Nelson Trujillo-Barreto %E%
0009 
0010 disp('Initialising volume');
0011 disp(' ');
0012 
0013 
0014 slice.X=X;
0015 slice.p=p;
0016 slice.k=size(X,2);
0017 
0018 T=size(X,1);
0019 slice.T=T;
0020 
0021 slice.XT=X';
0022 slice.XTX=slice.XT*X;
0023     
0024 for t=p+1:T,
0025     slice.dX(:,:,t-p) = X(t-1:-1:t-p,:);
0026 end
0027 
0028 [ux,dx,vx]=svd(X);
0029 ddx=diag(dx);
0030 svd_tol=max(ddx)*eps*slice.k;
0031 rank_X=sum(ddx > svd_tol);
0032 ddxm=diag(ones(rank_X,1)./ddx(1:rank_X));
0033 ddxm2=diag(ones(rank_X,1)./(ddx(1:rank_X).^2));
0034 slice.Xp=vx(:,1:rank_X)*ddxm*ux(:,1:rank_X)';
0035 slice.X2=vx(:,1:rank_X)*ddxm2*vx(:,1:rank_X)';
0036 
0037 k=slice.k;
0038 % Get lagged input covariances for updates
0039 slice.I.Gx=slice.X(p+1:T,:)'*slice.X(p+1:T,:);
0040 slice.I.xtx=slice.X(p+1:T,:)'*slice.X(p+1:T,:);
0041 for ki=1:k,
0042     for kj=1:k,
0043         dXki=squeeze(slice.dX(:,ki,:));
0044         dXkj=squeeze(slice.dX(:,kj,:));
0045         if slice.p==1
0046             % singleton dimension already transposed
0047             dXki=dXki';
0048             dXkj=dXkj';
0049         end
0050         Stmp = dXki*dXkj';
0051         slice.I.S((kj-1)*k+ki,:)=Stmp(:)';
0052         
0053         R1tmp= dXkj*slice.X(p+1:T,ki);
0054         slice.I.R1((kj-1)*k+ki,:)=R1tmp';
0055     end
0056 end
0057