lcd_trialfitls_holocalAR

 MAR parameter estimate with l2 reguralization for higher order local AR
 model (holocalAR model)

 * higher order local AR model 
  j(v,t) = sum_p a(v,p)*j(v,t-p) + sum_u b(v,u)*j(u,t-Delta(v,u)) + epsilon(t)
 
 [Input] 
 Zact :  current waveforms [Nv Nt Ntrial]
 Delta : delay matrix (sparse) with unit of time points 
         (i.e. all the elements must be integer representing dicrete time
         points)
   p    : AR order for local dynamics 
 
 [Output]
  A  : diagonal parts of MAR matrix (local dynamics)  [Nv,p]
  B  : non-digonal parts of MAR matrix (distant interactions) [Nv Nv]

 2017/09/23 OY

0001 function [A,B,Sig,AT,BT]=lcd_trialfitls_holocalAR(Zact,Delta,p,lambda)
0002 % MAR parameter estimate with l2 reguralization for higher order local AR
0003 % model (holocalAR model)
0004 %
0005 % * higher order local AR model
0006 %  j(v,t) = sum_p a(v,p)*j(v,t-p) + sum_u b(v,u)*j(u,t-Delta(v,u)) + epsilon(t)
0007 %
0008 % [Input]
0009 % Zact :  current waveforms [Nv Nt Ntrial]
0010 % Delta : delay matrix (sparse) with unit of time points
0011 %         (i.e. all the elements must be integer representing dicrete time
0012 %         points)
0013 %   p    : AR order for local dynamics
0014 %
0015 % [Output]
0016 %  A  : diagonal parts of MAR matrix (local dynamics)  [Nv,p]
0017 %  B  : non-digonal parts of MAR matrix (distant interactions) [Nv Nv]
0018 %
0019 % 2017/09/23 OY
0020 
0021 if nargin < 4
0022     lambda = 0;
0023 end
0024 
0025 [Nv,Nt,Ntr] = size(Zact);
0026 Dmax = max(Delta(:));
0027 Dmax = max([Dmax,p]);
0028 
0029 Zsum=sum(sum(abs(Zact),3),2);
0030 ixgood=find(Zsum ~= 0);
0031 Zgood = Zact(ixgood,:,:);
0032 Dgood = Delta(ixgood,ixgood);
0033 
0034 B = zeros(Nv,Nv);
0035 A = zeros(Nv,p);
0036 BT = zeros(Nv,Nv);
0037 AT = zeros(Nv,p);
0038 
0039 
0040 for vv = 1 : length(ixgood)
0041     
0042     fprintf('process %04d ...\n',ixgood(vv));
0043     
0044     
0045     dvec = round(Dgood(vv,:));
0046     iv = setdiff(find(dvec ~= 0), vv); % other than self loop
0047     
0048     % define X,y
0049     Ntt = Nt - Dmax;
0050     X = zeros(Ntt*Ntr, length(iv)+p);
0051     y = zeros(Ntt*Ntr,1);  
0052     for ntr = 1 : Ntr
0053         y((ntr-1)*Ntt+1:ntr*Ntt)=[Zgood(vv,Dmax+1:Nt,ntr)]';
0054         for aa = 1 : p
0055             X((ntr-1)*Ntt+1:ntr*Ntt,aa) = Zgood(vv,Dmax+1-aa:Nt-aa,ntr)';
0056         end
0057         nn = p+1;
0058         for uu = iv
0059             X((ntr-1)*Ntt+1:ntr*Ntt,nn) = Zgood(uu,Dmax+1-dvec(uu):Nt-dvec(uu),ntr)';
0060             nn = nn + 1;
0061         end
0062     end
0063     
0064    % least squares
0065    iXX = inv(X'*X + lambda*eye(length(iv)+p));
0066    Xy  = X'*y;
0067    
0068    w = iXX*Xy;
0069    e = y - X*w;
0070    sig = sum(e.^2)/(Ntt*Ntr - length(iv)-p);
0071    wT = w ./ sqrt(diag(iXX)) / sqrt(sig);
0072 
0073    if ixgood(vv) == 661, 
0074        display(vv)
0075    end
0076 
0077     % store reults
0078    A(ixgood(vv),:) = w(1:p);
0079    B(ixgood(vv),ixgood(iv)) = w(p+1:end);
0080    AT(ixgood(vv),:) = wT(1:p);
0081    BT(ixgood(vv),ixgood(iv)) = wT(p+1:end);
0082    Sig(ixgood(vv)) = sig;
0083 
0084     
0085     
0086     
0087     
0088     
0089 end
0090 
0091  
0092 
0093 
0094 
0095 
0096 
0097