lcd_fitl2reg_fukushima

 MAR parameter estimate with l2 reguralization for Fukushima model 

 * Original Fukushima's model 
  j(v,t) = a(v)*j(v,t-Delta(v,v)) + 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)
 lambda : regularization parameters 
          if set as a scalar, lambda is common to all the vertices
          if set [], optimal lambda is computed for each vertex  
          if set as a vector, A and B is computed for all the elements of
          lambda
 
 [Output]
  A  : diagonal parts of MAR matrix (local dynamics)  [Nv,1, length(lambda)]
  B  : non-digonal parts of MAR matrix (distant interactions) [Nv Nv length(lambda)]

 2016/05/02 O.Yamashita 
 * avoid increasing memory allocation
 2016/01/28 O.Yamashita

0001 function [A,B]=lcd_fitl2reg_fukushima(Zact,Delta,lambda)
0002 % MAR parameter estimate with l2 reguralization for Fukushima model
0003 %
0004 % * Original Fukushima's model
0005 %  j(v,t) = a(v)*j(v,t-Delta(v,v)) + sum_u b(v,u)*j(u,t-Delta(v,u)) + epsilon(t)
0006 %
0007 % [Input]
0008 % Zact :  current waveforms [Nv Nt]
0009 % Delta : delay matrix (sparse) with unit of time points
0010 %         (i.e. all the elements must be integer representing dicrete time
0011 %         points)
0012 % lambda : regularization parameters
0013 %          if set as a scalar, lambda is common to all the vertices
0014 %          if set [], optimal lambda is computed for each vertex
0015 %          if set as a vector, A and B is computed for all the elements of
0016 %          lambda
0017 %
0018 % [Output]
0019 %  A  : diagonal parts of MAR matrix (local dynamics)  [Nv,1, length(lambda)]
0020 %  B  : non-digonal parts of MAR matrix (distant interactions) [Nv Nv length(lambda)]
0021 %
0022 % 2016/05/02 O.Yamashita
0023 % * avoid increasing memory allocation
0024 % 2016/01/28 O.Yamashita
0025 
0026 
0027 [Nv,Nt] = size(Zact);
0028 Dmax = max(Delta(:));
0029 dself = Delta(1,1);  % delay of self loop
0030 
0031 if isempty(lambda) | length(lambda) == 1
0032     B = zeros(Nv,Nv);
0033     A = zeros(Nv,1);
0034 else
0035     B = zeros(Nv,Nv,length(lambda));
0036     A = zeros(Nv,1,length(lambda));
0037 end
0038 
0039 for vv = 1 : Nv
0040     fprintf('process %04d ...\n',vv);
0041     
0042     y = [Zact(vv,Dmax+1:Nt)]';
0043        
0044     dvec = Delta(vv,:);
0045     iv = setdiff(find(dvec ~= 0), vv);
0046     
0047     X = zeros(Nt-Dmax, length(iv)+1);
0048     X(:,1) = Zact(vv,Dmax+1-dself:Nt-dself)';
0049     
0050     nn = 2;
0051     for uu = iv
0052         X(:,nn) = Zact(uu,Dmax+1-dvec(uu):Nt-dvec(uu))';
0053         nn = nn + 1;
0054     end
0055        
0056     
0057     if isempty(lambda) | length(lambda) == 1
0058         [lam,sig,val,a] = wmn_solution(y,X,[],[],struct('flag_plot',0,'lambda',lambda));
0059         A(vv,:) = a(1);    
0060         B(vv,iv) = a(2:end);
0061     else
0062         [aa] = wmn_solution_fixreg(y, X, [], [], lambda, struct('flag_plot',0));
0063         A(vv,1,:) = aa(1,:);
0064         B(vv,iv,:) = aa(2:end,:);
0065     end
0066     
0067 end
0068 
0069 %
0070 % figure,
0071 % subplot(2,1,1)
0072 % plot(A);
0073 % subplot(2,1,2)
0074 % plot(B');
0075 %
0076 %
0077 % ixeff = find(abs(mean(A,2)) > 0.05) ;
0078 % figure
0079 % subplot(2,1,1)
0080 % plot(A(ixeff,:)),
0081 % set(gca, 'xtick', [1:length(ixeff)], 'xticklabel', num2str(ixeff));
0082 % subplot(2,1,2)
0083 % imagesc(B(ixeff,ixeff),[-0.05,0.05]);
0084 % set(gca, 'xtick', [1:length(ixeff)], 'xticklabel', num2str(ixeff), 'ytick', [1:length(ixeff)], 'yticklabel', num2str(ixeff));
0085 % cbar
0086 
0087 
0088 
0089 
0090 
0091