lcd_fitrbfmar_fukushima

 [Input] 
  X :
 Delta :
  ti : center time points of rbf kernel
 s2  : rbf kernel width 

 2016/03/22 O.Yamashita

0001 function [A,B,C,Xc]= lcd_fitrbfmar_fukushima(X,Delta,ti,s2, varargin)
0002 %
0003 % [Input]
0004 %  X :
0005 % Delta :
0006 %  ti : center time points of rbf kernel
0007 % s2  : rbf kernel width
0008 %
0009 % 2016/03/22 O.Yamashita
0010 
0011 
0012 parm = finputcheck(varargin, ...
0013     {'nlearn'     ,'integer',  [1 inf],  100;...
0014      'lam0'       ,'real'   ,  []     ,  1e-4;...
0015      'beta0'      ,'real'   ,  [0 inf],  1;...
0016      'd0'         ,'real'   ,  [0 inf],  1;...   
0017      });
0018  
0019 if ~isstruct(parm)
0020    error(parm);
0021 end
0022 
0023 Nlearn = parm.nlearn;
0024 lam0  = parm.lam0;
0025 beta0 = parm.beta0;
0026 d0    = parm.d0;
0027 
0028 dmax = max(Delta(:));
0029 [Nv,Nt]=size(X);
0030 
0031 % make input regressors
0032 for d = 1 : length(ti)
0033     Xc(:,d) = exp(-([dmax+1:Nt]'-ti(d)).^2/s2);
0034 end
0035 XcXc=Xc'*Xc;
0036 
0037 
0038 A = zeros(Nv,1);
0039 B = zeros(Nv,Nv);
0040 C = zeros(Nv,length(ti));
0041 
0042 for vv = 1 : Nv
0043 
0044     y = X(vv,dmax+1:Nt)';
0045     ix = find(Delta(vv,:)~=0);
0046 
0047     
0048     Xb = [];
0049     for uu = 1 : length(ix)
0050         if ix(uu) == vv
0051             Xa = X(vv,dmax+1-Delta(vv,vv):Nt-Delta(vv,vv))';
0052         else
0053             Xb = [Xb X(ix(uu),dmax+1-Delta(vv,ix(uu)):Nt-Delta(vv,ix(uu)))'];
0054         end
0055     end
0056 
0057     if isempty(Xb)
0058         Xb = 0;
0059         lam0 = 1;
0060         emptyXb = 1;
0061     else
0062         emptyXb = 0;
0063     end
0064     XbXb = Xb'*Xb;
0065     XaXa = Xa'*Xa;
0066  
0067 
0068     %
0069     % learning
0070     %
0071 
0072     lam = lam0;
0073     beta = beta0;
0074     a    = zeros(size(Xa,2),1);
0075     c    = zeros(size(Xc,2),1);
0076     d    = d0*ones(size(Xc,2),1);
0077 
0078     for nn = 1 : Nlearn
0079 
0080 
0081         Pb = beta*XbXb+lam*eye(size(Xb,2));  % Pb = beta*Xb'*Xb+lam*eye(size(Xb,2));
0082         Sb = inv(Pb);
0083         b = Sb *(Xb'*(y-Xa*a-Xc*c)*beta);
0084 
0085         Pa = beta*XaXa; %Pa = beta*Xa'*Xa;
0086         Sa = inv(Pa);
0087         a  = Sa*(Xa'*(y-Xb*b-Xc*c)*beta);
0088 
0089         Pc = beta*XcXc+diag(d); % Pc = beta*Xc'*Xc+diag(d);
0090         Sc = inv(Pc);
0091         c = Sc*(Xc'*(y-Xa*a-Xb*b)*beta);
0092 
0093 %        d = 1./(diag(Sc)+c.^2);
0094         d = (1 - c.^2)./diag(Sc);
0095         lam = length(b)/(b'*b+trace(Sb));
0096 %        beta = (Nt-dmax)/(sum((y-Xa*a-Xb*b-Xc*c).^2)+trace(Xa'*Xa*Pa)+trace(Xb'*Xb*Pb)+trace(Xc'*Xc*Pc));
0097 %        beta = (Nt-dmax)/(sum((y-Xa*a-Xb*b-Xc*c).^2)+trace(XaXa*Pa)+trace(XbXb*Pb)+trace(XcXc*Pc));
0098         beta = (Nt-dmax)/(sum((y-Xa*a-Xb*b-Xc*c).^2)+sum(sum(XaXa.*Pa'))+sum(sum(XbXb.*Pb'))+sum(sum(XcXc.*Pc')));
0099 
0100     end
0101 
0102     
0103     if emptyXb
0104         A(vv) = a;
0105         C(vv,:) = c;
0106     else
0107         A(vv) = a;
0108         B(vv,setdiff(ix,vv)) = b;
0109         C(vv,:) = c;
0110     end
0111 
0112 
0113 end