nt_tsregress

[z,idx]=nt_tsregress(x,y,shifts,xw,yw,keep,threshold) - time-shift regression

  z: part of x modeled by time-shifted y
  idx: x(idx) maps to z

  x: data to model
  y: regressor
  shifts: array of shifts to apply (default: [0])
  xw: weights to apply to x
  yw: weights to apply to y
  keep: number of components of shifted regressor PCs to keep (default: all)
  threshold: discard PCs with eigenvalues below this (default: 0)

 Data X are regressed on time-shifted versions of Y. X and Y are initially 
 time-aligned, but because of the shifts, Z is shorter than X.  Z is
 time-aligned with X(IDX).

0001 function [z,idx]=nt_tsregress(x,y,shifts,xw,yw,keep,threshold)
0002 %[z,idx]=nt_tsregress(x,y,shifts,xw,yw,keep,threshold) - time-shift regression
0003 %
0004 %  z: part of x modeled by time-shifted y
0005 %  idx: x(idx) maps to z
0006 %
0007 %  x: data to model
0008 %  y: regressor
0009 %  shifts: array of shifts to apply (default: [0])
0010 %  xw: weights to apply to x
0011 %  yw: weights to apply to y
0012 %  keep: number of components of shifted regressor PCs to keep (default: all)
0013 %  threshold: discard PCs with eigenvalues below this (default: 0)
0014 %
0015 % Data X are regressed on time-shifted versions of Y. X and Y are initially
0016 % time-aligned, but because of the shifts, Z is shorter than X.  Z is
0017 % time-aligned with X(IDX).
0018 
0019 if nargin<2; error('!'); end
0020 if nargin<3||isempty(shifts); shifts=[0]; end
0021 if nargin<4; xw=[]; end
0022 if nargin<5; yw=[]; end
0023 if nargin<6; keep=[]; end
0024 if nargin<7; threshold=[]; end
0025 
0026 if size(x,1) ~= size(y,1); error('!'); end
0027 
0028 % shifts must be non-negative
0029 mn=min(shifts);
0030 if mn<0; 
0031     shifts=shifts-mn; 
0032     x=x(-mn+1:end,:,:);
0033     y=y(-mn+1:end,:,:);
0034 end
0035 nshifts=numel(shifts);
0036 
0037 % % flag outliers in x and y
0038 % if ~isempty(toobig1) || ~isempty(toobig2)
0039 %     xw=nt_find_outliers(x,toobig1,toobig2);
0040 %     yw=nt_find_outliers(y,toobig1,toobig2);
0041 % else
0042 %     xw=[];yw=[];
0043 %     %xw=ones(size(x)); yw=ones(size(y));
0044 % end
0045 
0046 % subtract weighted means
0047 
0048 if ndims(x)==3    
0049     [Mx,Nx,Ox]=size(x);
0050     [My,Ny,Oy]=size(y);
0051     x=nt_unfold(x);
0052     y=nt_unfold(y);
0053     [x,xmn]=nt_demean(x,xw);
0054     [y,ymn]=nt_demean(y,yw);
0055     x=nt_fold(x,Mx);
0056     y=nt_fold(y,My);
0057 else
0058     [x,xmn]=nt_demean(x,xw);
0059     [y,ymn]=nt_demean(y,yw);
0060 end
0061 
0062 
0063 % covariance of y
0064 [cyy,totalweight]=nt_cov(y,shifts',yw);
0065 cyy=cyy./totalweight;
0066 
0067 % cross-covariance of x and y
0068 [cxy, totalweight]=nt_cov2(x,y,shifts',xw,yw);
0069 disp('!!!!!!!!!   WARNING: calling obsolete code  !!!!!!!!!!!!!!!!');
0070 %[cxy, totalweight]=nt_xcov(x,y,shifts',xw,yw);
0071 cxy=cxy./totalweight;
0072 
0073 % regression matrix
0074 r=nt_regcov(cxy,cyy,keep,threshold);
0075     
0076 % regression
0077 if ndims(x)==3
0078     x=nt_unfold(x);
0079     y=nt_unfold(y);
0080     [m,n]=size(x);
0081     mm=m-max(shifts);
0082     z=zeros(size(x));
0083     for k=1:nshifts
0084         kk=shifts(k);
0085         idx1=kk+1:kk+mm;
0086         idx2=k+(0:size(y,2)-1)*nshifts;
0087         z(1:mm,:)=z(1:mm,:)+y(idx1,:)*r(idx2,:);
0088     end
0089     z=nt_fold(z,Mx);
0090     z=z(1:end-max(shifts),:,:);
0091 else
0092     [m,n]=size(x);
0093     z=zeros(m-max(shifts),n);
0094     for k=1:nshifts
0095         kk=shifts(k);
0096         idx1=kk+1:kk+size(z,1);
0097         %idx2=k*size(y,2)+1:(k+1)*size(y,2);
0098         idx2=k+(0:size(y,2)-1)*nshifts;
0099         z=z+y(idx1,:)*r(idx2,:);
0100     end
0101 end
0102 
0103 % idx allows x to be aligned with z
0104 offset=max(0,-mn);
0105 idx=offset+1:offset+size(z,1);
0106 
0107 %elseif ndims(x)==3
0108     
0109 %end
0110