home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / ARMAX.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  6.3 KB  |  190 lines
  1. function th=armax(z,nn,maxiter,tol,lim,maxsize,Tcap)
  2. %ARMAX  Computes the prediction error estimate of an ARMAX model.
  3. %
  4. %    TH = armax(Z,NN)
  5. %
  6. %    TH: returned as  the estimated parameters of the ARMAX model
  7. %    A(q) y(t) = B(q) u(t-nk) + C(q) e(t)
  8. %    along with estimated covariances and structure information.
  9. %    For the exact format of TH, see HELP THETA.
  10. %
  11. %    Z : The output-input data Z=[y u], with y and u being column vectors.
  12. %    For a time-series, Z=y only. The routine does not work for multi-
  13. %    input systems. Use PEM for that case.
  14. %
  15. %    NN: Initial value and structure information. When no initial parameter
  16. %    estimates are available enter NN as
  17. %    NN=[na nb nc nk], the orders and delay of the above model,or as
  18. %    NN=[na nc] for the time-series case (ARMA-model). With an initial
  19. %    estimate available in THI, a theta matrix of standard format, enter
  20. %    NN=THI. Then the criterion minimization is initialized at THI.
  21. %
  22. %    Some parameters associated with the algorithm are accessed by
  23. %    TH = armax(Z,NN,maxiter,tol,lim,maxsize,T)
  24. %    See HELP AUXVAR for an explanation of these and their default values.
  25.  
  26. %    L. Ljung 10-1-86,12-09-91
  27. %    Copyright (c) 1986-92 by the MathWorks, Inc.
  28. %    All Rights Reserved.
  29.  
  30. %
  31. % *** Set up default values ***
  32. [Ncap,nz]=size(z);
  33. [nr,cc]=size(nn);   nu=nz-1;
  34. maxsdef=idmsize(Ncap);
  35.  
  36. if nargin<7, Tcap=[];end
  37. if nargin<6, maxsize=[];end
  38. if nargin<5, lim=[];end
  39. if nargin<4, tol=[];, end
  40. if nargin<3, maxiter=[]; end
  41. if isempty(Tcap),if nr>1,Tcap=gett(nn);else Tcap=1;end,end
  42. if isempty(maxsize),maxsize=maxsdef;end
  43. if isempty(lim),lim=1.6;end,if isempty(tol),tol=0.01;end
  44. if isempty(maxiter),maxiter=10;end
  45. if Tcap<0,Tcap=1;end, if maxsize<0,maxsize=maxsdef;end, if lim<0,lim=1.6;end
  46. if tol<0,tol=0.01;end,if maxiter<0, maxiter=10;end
  47.  
  48. % *** Do some consistency tests ***
  49.  
  50.  
  51. if nz>Ncap, error('The data should be organized in column vectors')
  52. return,end
  53. if nu>1, error('This routine only works for single-input systems. Use PEM instead!')
  54. return,end
  55. if nr==1 & cc~=2*nu+2, 
  56.        disp('Incorrect number of orders specified:')
  57.        disp('For a time series  nn=[na nc]')
  58.        disp('For a SISO system  nn=[na nb nc nk]')
  59.        error('see above')
  60. return,end
  61.  
  62. % *** if nn has more than 1 row, it is a theta matrix
  63. %      and we jump directly to the minimization ***
  64. %
  65. if nr>1, nu=nn(1,3);
  66.    if nu>1, error('This routine only works for single input systems. Use PEM instead')
  67.    return,end
  68.    if nu==1,  na=nn(1,4);nb=nn(1,5);nc=nn(1,6);nk=nn(1,9);
  69.    else na=nn(1,4);nb=0;nc=nn(1,5);nk=0;end
  70.          n=na+nb+nc; t=nn(3,1:n).';
  71.          if nc>0, c=[1 t(na+nb+1:n).'];else c=1;,end
  72.          if na>0,a=[1 t(1:na).'];else a=1;,end
  73.          if nb>0,b=[zeros(1,nk) t(na+1:na+nb).'];,else b=0;,end
  74.          ni=max([na nb+nk-1 nc]);
  75.      e=pefilt(a,c,z(:,1)); if nu==1, e=e-pefilt(b,c,z(:,2),e(1:ni));end
  76. end
  77. if nr==1,
  78.    if nu==0,
  79.           na=nn(1); nc=nn(2); nb=0; nk=0; n=na+nc;
  80.           v=z; narx=na; nd=na;
  81.    end
  82.    if nu==1,
  83.           na=nn(1); nb=nn(2); nc=nn(3); nk=nn(4); n=sum(nn(1:3));
  84.           narx=[na nb nk]; nd=0;
  85.    end
  86.    ni=max([na nb+nk-1 nc]);
  87.    t=arx(z,narx,maxsize,Tcap,0);
  88.    if na>0, a=[1 t(1:na).'];else a=1;end
  89.    if nb>0, b=[zeros(1,nk) t(na+1:na+nb).'];end
  90.    if nc==0  % then we use the arx-estimate for initial condition
  91.        e=pefilt(a,[1],z(:,1));
  92.        if nu==1, e=e-pefilt(b,[1],z(:,2),e(1:ni));end
  93.        c=1;
  94.    end
  95.  
  96.    if nc>0 %then we compute the IV-estimate in case nu==1
  97.       if nu==1
  98.          a=fstab(a);
  99.          t=iv(z,narx,a,b,maxsize,Tcap,0);
  100.          if na>0,a=[1 t(1:na).'];else a=1;end
  101.          b=[zeros(1,nk) t(na+1:na+nb).'];
  102.          v=filter(a,1,z(:,1))-filter(b,1,z(:,2));v(1:ni)=zeros(ni,1);
  103.       end
  104. %
  105. %  *** Determine the initial C-estimate by first building a
  106. %      high order AR-model of the residuals v and the use the 
  107. %      LS-method with the new residuals as inputs ***
  108. %          
  109.         ord=min(na+nb+4*nc,length(v)/3);
  110.     t1=arx(v,ord,maxsize,Tcap,0);
  111.         e=pefilt([1 t1.'],1,v);
  112.         t2=arx([v e],[nd nc 1],maxsize,Tcap,0);
  113. %
  114. %         ** test the stability of the C-estimate **
  115. %
  116.         c=[1 t2(nd+1:nd+nc).']; c=fstab(c);
  117.         t2(nd+1:nd+nc)=c(2:nc+1).';
  118.    
  119.      if nu==1, t=[t;t2]; else t=t2;end
  120. %
  121.      if nd==0, d=1; else d=[1 t2(1:nd).'];end
  122.      e=pefilt(d,c,v,zeros(1,ni));
  123.   end
  124. end
  125. % *** Display initial estimate ***
  126. %
  127. V=e'*e/(Ncap-ni);
  128. clc, disp(['  INITIAL ESTIMATE'])
  129. disp(['Current loss:' num2str(V)])
  130. disp(['th-vector'])
  131. disp(t)
  132. %
  133. % *** start minimizing ***
  134. %
  135. % ** determine limit for robustification **
  136.  if lim~=0, lim=median(abs(e-median(e)))*lim/0.7;end
  137.  if lim==0,el=e;else
  138.  [ne,me]=size(e);
  139.  ll=ones(ne,me)*lim;la=abs(e)+eps*ll;el=e.*(min(la,ll)./la); clear ll,clear la,end
  140.  g=ones(n,1); l=0; st=0; nmax=max([na nb+nk-1 nc]);
  141. %
  142. % ** the minimization loop **
  143. %
  144. while [norm(g)>tol l<maxiter st~=1]
  145.       l=l+1;
  146. %     * compute gradient *
  147.       yf=filter(-1,c,z(:,1)); ef=filter(1,c,e);
  148.       if nu==1, uf=filter(1,c,z(:,2));end
  149.   M=floor(maxsize/n);
  150.   R=zeros(n);F=zeros(n,1);
  151.   for k1=nmax:M:Ncap-1
  152.       jj=(k1+1:min(Ncap,k1+M));
  153.       psi=zeros(length(jj),n);
  154.       for k=1:na, psi(:,k)=yf(jj-k);end
  155.       for k=1:nb, psi(:,na+k)=uf(jj-k-nk+1);end
  156.       for k=1:nc, psi(:,k+nb+na)=ef(jj-k);end
  157.       R=R+psi'*psi; F=F+psi'*el(jj);
  158.   end
  159.    if Ncap>M, g=R\F; else g=psi\el(jj);end,grad=F;
  160. %
  161. %     * search along the g-direction *
  162. %
  163.       [t1,e,el,V1,c,st]=searchax(z,t,g,lim,V,na,nb,nc,nk,ni);
  164. if st==1,      [t1,e,el,V1,c,st]=searchax(z,t,grad/trace(R)*length(R),lim,V,na,nb,nc,nk,ni);end
  165.       home
  166.       disp(['      ITERATION # ' int2str(l)])
  167.       disp(['Current loss: ' num2str(V1) '  Previous loss: ' num2str(V)])
  168.       disp(['Current th  prev. th   gn-dir'])
  169.       disp([t1 t g])
  170.       disp(['Norm of gn-vector: '  num2str(norm(g))])
  171.  
  172.       if st==1, disp(['No improvement of the criterion possible along the search direction'])
  173.          disp(['Iterations therefore terminated']),end
  174.       t=t1; V=V1;
  175. end
  176. th=zeros(3+n,max([6+3*nu 7 n]));
  177. V=e'*e/(length(e)-ni);
  178. if nu==1
  179.    th(1,1:9)=[V Tcap 1 na nb nc 0 0 nk];
  180.    else th(1,1:6)=[V Tcap 0 na nc 0];end
  181. th(2,1)=V*(1+n/Ncap)/(1-n/Ncap);
  182. ti=fix(clock); ti(1)=ti(1)/100;
  183. th(2,2:6)=ti(1:5);
  184. th(2,7)=7;
  185. th(3,1:n)=t.';
  186. if maxiter==0,return,end
  187. if Ncap>M, PP=inv(R); else PP=inv(psi'*psi);end
  188. th(4:3+n,1:n)=V*PP;
  189.  
  190.