lcd_fitl2reg2step_holocalAR

 Two-step 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]
 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 
 lambda : a regularization parameter common to all the vertices
          if set [], optimal lambda is computed for each vertex  
 
 [Output]
  A  : diagonal parts of MAR matrix (local dynamics)  [Nv,p]
  B  : non-digonal parts of MAR matrix (distant interactions) [Nv Nv]

 2016/02/01 O.Yamashita

0001 function [A,B,Lam,Sig]=lcd_fitl2reg2step_holocalAR(Zact,Delta,p,lambda)
0002 % Two-step 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]
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 % lambda : a regularization parameter common to all the vertices
0015 %          if set [], optimal lambda is computed for each vertex
0016 %
0017 % [Output]
0018 %  A  : diagonal parts of MAR matrix (local dynamics)  [Nv,p]
0019 %  B  : non-digonal parts of MAR matrix (distant interactions) [Nv Nv]
0020 %
0021 % 2016/02/01 O.Yamashita
0022 
0023 
0024 [Nv,Nt] = size(Zact);
0025 Dmax = max(Delta(:));
0026 Dmax = max([Dmax,p]);
0027 
0028 B = zeros(Nv,Nv);
0029 A = zeros(Nv,p);
0030 
0031 
0032 for vv = 1 : Nv
0033     
0034     fprintf('process %04d ...\n',vv);
0035     
0036     y = [Zact(vv,Dmax+1:Nt)]';
0037     
0038     %
0039     % distant regressors
0040     %
0041     
0042     X = []; 
0043     dvec = round(Delta(vv,:));
0044     iv = setdiff(find(dvec ~= 0), vv); % other than self loop
0045     
0046     for uu = iv
0047         X = [X, Zact(uu,Dmax+1-dvec(uu):Nt-dvec(uu))'];
0048     end
0049     
0050     if ~isempty(X)
0051         [lam,sig,val,a] = wmn_solution(y,X,[],[],struct('flag_plot',0,'lambda',lambda));
0052         B(vv,iv) = a;
0053         e = y - X*a;  % residuals
0054     else
0055         e = y;
0056     end
0057 
0058     %
0059     % local regressors
0060     %
0061     X = [];
0062     for aa = 1 : p
0063     X = [X Zact(vv,Dmax-aa+1:Nt-aa)'];
0064     end
0065     [lam,sig,val,a] = wmn_solution(e,X,[],[],struct('flag_plot',0,'lambda',lambda));
0066 
0067     A(vv,:) = a;    
0068 
0069     
0070     Lam(vv) = lam;
0071     Sig(vv) = sig;
0072     
0073 end
0074 
0075 
0076 % figure,
0077 % subplot(2,1,1)
0078 % plot(A);
0079 % subplot(2,1,2)
0080 % plot(B');
0081 %
0082 
0083 
0084 
0085 
0086 
0087