home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / RBJ.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.8 KB  |  98 lines
  1. function [thm,yhat,p,phi,psi] = rbj(z,nn,adm,adg,th0,p0,phi,psi)
  2. %RBJ    Computes estimates recursively for a BOX-JENKINS model using the
  3. %    Recursive Prediction Error Method
  4. %
  5. %    [THM,YHAT] = rbj(Z,NN,adm,adg)
  6. %
  7. %    Z: The output-input data z=[y u] (single input only!)
  8. %    NN : NN=[nb nc nd nf nk], The orders and delay of a general
  9. %         input-output model (see HELP BJ)
  10. %    adm: Adaptation mechanism. adg: Adaptation gain
  11. %     adm='ff', adg=lam:  Forgetting factor algorithm, with forg factor lam
  12. %     adm='kf', adg=R1: The Kalman filter algorithm with R1 as covariance 
  13. %        matrix of the parameter changes per time step
  14. %     adm='ng', adg=gam: A normalized gradient algorithm, with gain gam
  15. %     adm='ug', adg=gam: An Unnormalized gradient algorithm with gain gam
  16. %    THM: The resulting estimates. Row k contains the estimates "in alpha-
  17. %         betic order" corresponding to data up to time k (row k in Z)
  18. %    YHAT: The predicted values of the outputs. Row k corresponds to time k
  19. %    Initial value of parameters(TH0) and of "P-matrix" (P0) can be given by
  20. %    [THM,YHAT,P] = rbj(Z,NN,adm,adg,TH0,P0)
  21. %    Initial and last values of auxiliary data vectors phi and psi are 
  22. %    obtained by [THM,YHAT,P,phi,psi]=rbj(Z,NN,adm,adg,TH0,P0,phi0,psi0). 
  23. %    
  24. %    See also RARX, RARMAX, ROE, RPEM and RPLR.
  25.  
  26. %    L. Ljung 10-1-89
  27. %    Copyright (c) 1989-90 by the MathWorks, Inc.
  28. %    All Rights Reserved.
  29.  
  30. [nz,ns]=size(z);
  31. if ns<=1,error('This routine requires an input. For a time series, use RARMAX or RARX instead!'),end
  32. if ns>2,error('This  routine is for single input only. Use RPEM instead!'),end
  33. if length(nn)~=5,error('Incorrect number of orders specified! nn = [nb nc nd nf nk]'),end
  34. nb=nn(1);nc=nn(2);nd=nn(3);nf=nn(4);nk=nn(5);nu=1;
  35. if nk<1,error('Sorry, this routine requires nk>0; Shift input sequence if necessary!'),end
  36. d=sum(nn(1:4));
  37. if ns>2,error('Sorry, this routine is for single input only!'),end
  38. ng=nf+nc;
  39. nbm=max([nb+nk-1,ng,nd]);ndm=max(nd,nc);nfm=max([nf,ng,nd]);
  40. tic=nb+1:nb+nc;tif=nb+nc+nd+1:d;
  41. tib=1:nb;tid=nb+nc+1:nb+nc+nd;
  42. ib=nk:nb+nk-1;ibg=1:ng;ibd=nk:nk+nd-1;
  43. ic=nbm+1:nbm+nc;
  44. id=nbm+nc+1:nbm+nc+nd;idc=nbm+nc+1:nbm+nc+nc;
  45. iff=nbm+nc+ndm+1:nbm+nc+ndm+nf;ifg=nbm+nc+ndm+1:nbm+nc+ndm+ng;
  46. ifd=nbm+nc+ndm+1:nbm+nc+ndm+nd;
  47. dm=nfm+nbm+nc+ndm;
  48. i=[ib ic id iff];
  49. if nargin<8, psi=zeros(dm,1);end
  50. if nargin<7, phi=zeros(dm,1);end
  51. if nargin<6, p0=10000*eye(d);end
  52. if nargin<5, th0=eps*ones(d,1);end
  53. if isempty(psi),psi=zeros(dm,1);end
  54. if isempty(phi),phi=zeros(dm,1);end
  55. if isempty(p0),p0=10000*eye(d);end
  56. if isempty(th0),th0=eps*ones(d,1);end
  57. if length(th0)~=d, error('The length of th0 must equal the number of estimated parameters!'),end
  58. [th0nr,th0nc]=size(th0);if th0nr<th0nc, th0=th0';end
  59. p=p0;th=th0;
  60. if adm(1)=='f', R1=zeros(d,d);lam=adg;end
  61. if adm(1)=='k', [sR1,SR1]=size(adg);
  62.      if sR1~=d | SR1~=d,error('The R1 matrix should be a square matrix with dimension equal to number of parameters!'),end
  63.      R1=adg;lam=1;
  64. end
  65. if adm(2)=='g', grad=1;else grad=0;end
  66.  
  67. for kcou=1:nz
  68. yh=phi(i)'*th;
  69. epsi=z(kcou,1)-yh;
  70. if ~grad,K=p*psi(i)/(lam + psi(i)'*p*psi(i));
  71.          p=(p-K*psi(i)'*p)/lam+R1;
  72. else K=adg*psi(i);end
  73. if adm(1)=='n', K=K/(eps+psi(i)'*psi(i));end
  74. th=th+K*epsi;
  75. c=fstab([1;th(tic)])';f=fstab([1;th(tif)])';d=[1;th(tid)];
  76. th(tic)=c(2:nc+1);th(tif)=f(2:nf+1);g=conv(f,c);%HIT*********
  77. w=th([tib tif])'*phi([ib iff]);
  78. util=d'*[z(kcou,2);phi(ibd)]-g'*[0;psi(ibg)];
  79. if nf>0
  80. wtil=d'*[w;-phi(ifd)]+g'*[0;psi(ifg)];
  81. end
  82. v=z(kcou,1)-w;
  83. epsilon=v-th([tic tid])'*phi([ic id]);
  84.  
  85. if nc>0
  86. epstil=c'*[epsilon;-psi(ic)];end
  87. if nd>0
  88. vtil=c'*[v;psi(idc)];
  89. end
  90. phi(2:dm)=phi(1:dm-1);psi(2:dm)=psi(1:dm-1);
  91.  
  92. if nb>0,phi(1)=z(kcou,2);psi(1)=util;end 
  93. if nc>0,phi(ic(1))=epsilon;psi(ic(1))=epstil;end
  94. if nd>0,phi(id(1))=-v;psi(id(1))=-vtil;end
  95. if nf>0,phi(iff(1))=-w;psi(iff(1))=-wtil;end
  96. thm(kcou,:)=th';yhat(kcou)=yh;
  97. end
  98. yhat=yhat';