home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / RARMAX.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.6 KB  |  93 lines
  1. function [thm,yhat,p,phi,psi] = rarmax(z,nn,adm,adg,th0,p0,phi,psi)
  2. %RARMAX    Computes estimates recursively for an ARMAX model using the
  3. %    Recursive Prediction Error Method
  4. %
  5. %    [THM,YHAT] = rarmax(Z,NN,adm,adg)
  6. %
  7. %    Z: The output-input data z=[y u] (single input only!)
  8. %    NN : NN=[na nb nc nk], The orders and delay of an ARMAX
  9. %         input-output model (see HELP ARMAX)
  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 output. 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] = rarmax(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]=rarmax(Z,NN,adm,adg,TH0,P0,phi0,psi0). 
  23. %    
  24. %    See also RARX, ROE, RBJ, RPEM and RPLR.
  25.  
  26. %    L. Ljung 10-1-89
  27. %    Copyright (c) 1989 by the MathWorks, Inc.
  28. %    All Rights Reserved.
  29.  
  30.  
  31. [nz,ns]=size(z);[ordnr,ordnc]=size(nn);
  32. if ns>2,error('This  routine is for single input only. Use RPEM instead!'),end
  33. if ns==1, if ordnc~=2;error('For a time series nn should be [na nc]!'),end
  34. else if ordnc~=4, error('the argument nn should be [na nb nc nk]!'),end,end
  35. if ns==1,na=nn(1);nb=0;nc=nn(2);nk=1;
  36. else na=nn(1);nb=nn(2);nc=nn(3);nk=nn(4);nu=1;
  37. end
  38. if nk<1,error('Sorry, this routine requires nk>0; Shift input sequence if necessary!'),end
  39. d=na+nb+nc;
  40. if ns>2,error('Sorry, this routine is for single input only!'),end
  41. if ns==1,nb=0;end
  42. if nb==0,nk=1;end
  43. nam=max([na,nc]);nbm=max([nb+nk-1,nc]);
  44. tic=na+nb+1:na+nb+nc;
  45. ia=1:na;iac=1:nc;
  46. ib=nam+nk:nam+nb+nk-1;ibc=nam+1:nam+nc;
  47. ic=nam+nbm+1:nam+nbm+nc;
  48.  
  49. iia=1:nam-1;iib=nam+1:nam+nbm-1;iic=nam+nbm+1:nam+nbm+nc-1;
  50. dm=nam+nbm+nc; if nb==0,iib=[];end
  51. ii=[iia iib iic];i=[ia ib ic];
  52.  
  53. if nargin<8, psi=zeros(dm,1);end
  54. if nargin<7, phi=zeros(dm,1);end
  55. if nargin<6, p0=10000*eye(d);end
  56. if nargin<5, th0=eps*ones(d,1);end
  57. if isempty(psi),psi=zeros(dm,1);end
  58. if isempty(phi),phi=zeros(dm,1);end
  59. if isempty(p0),p0=10000*eye(d);end
  60. if isempty(th0),th0=eps*ones(d,1);end
  61. if length(th0)~=d, error('The length of th0 must equal the number of estimated parameters!'),end
  62. [th0nr,th0nc]=size(th0);if th0nr<th0nc, th0=th0';end
  63. p=p0;th=th0;
  64. if adm(1)=='f', R1=zeros(d,d);lam=adg;end
  65. if adm(1)=='k', [sR1,SR1]=size(adg);
  66.      if sR1~=d | SR1~=d,error('The R1 matrix should be a square matrix with dimension equal to number of parameters!'),end
  67.      R1=adg;lam=1;
  68. end
  69. if adm(2)=='g', grad=1;else grad=0;end
  70.  
  71. for kcou=1:nz
  72. yh=phi(i)'*th;
  73. epsi=z(kcou,1)-yh;
  74. if ~grad,K=p*psi(i)/(lam + psi(i)'*p*psi(i));
  75.          p=(p-K*psi(i)'*p)/lam+R1;
  76. else K=adg*psi(i);end
  77. if adm(1)=='n', K=K/(eps+psi(i)'*psi(i));end
  78. th=th+K*epsi;
  79. if nc>0,c=fstab([1;th(tic)])';else c=1;end
  80. th(tic)=c(2:nc+1);
  81. epsilon=z(kcou,1)-phi(i)'*th;
  82. if nb>0,zb=[z(kcou,2),-psi(ibc)'];else zb=[];end
  83. ztil=[[z(kcou,1),psi(iac)'];zb;[epsilon,-psi(ic)']]*c;
  84.  
  85. phi(ii+1)=phi(ii);psi(ii+1)=psi(ii);
  86. if na>0,phi(1)=-z(kcou,1);psi(1)=-ztil(1);end
  87. if nb>0,phi(nam+1)=z(kcou,2);psi(nam+1)=ztil(2);end 
  88. if nb==0,zc=ztil(2);else zc=ztil(3);end
  89. if nc>0,phi(nam+nbm+1)=epsilon;psi(nam+nbm+1)=zc;end
  90.  
  91. thm(kcou,:)=th';yhat(kcou)=yh;
  92. end
  93. yhat=yhat';