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