home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / RARX.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.8 KB  |  77 lines
  1. function [thm,yhat,p,phi] = rarx(z,nn,adm,adg,th0,p0,phi)
  2. %RARX    Computes estimates recursively for an ARX model using the
  3. %    Recursive Least Squares Method
  4. %
  5. %    [THM,YHAT] = rarx(Z,NN,adm,adg)
  6. %
  7. %    Z: The output-input data z=[y u]
  8. %    NN : NN=[na nb nk], The orders and delay of an ARX model (see HELP ARX)
  9. %    adm: Adaptation mechanism. adg: Adaptation gain
  10. %     adm='ff', adg=lam:  Forgetting factor algorithm, with forg factor lam
  11. %     adm='kf', adg=R1: The Kalman filter algorithm with R1 as covariance 
  12. %        matrix of the parameter changes per time step
  13. %     adm='ng', adg=gam: A normalized gradient algorithm, with gain gam
  14. %     adm='ug', adg=gam: An Unnormalized gradient algorithm with gain gam
  15. %    THM: The resulting estimates. Row k contains the estimates "in alpha-
  16. %         betic order" corresponding to data up to time k (row k in Z)
  17. %    YHAT: The predicted values of the outputs. Row k corresponds to time k.
  18. %    Initial value of parameters(TH0) and of "P-matrix" (P0) can be given by
  19. %    [THM,YHAT,P] = rarx(Z,NN,adm,adg,TH0,P0)
  20. %    Initial and last values of auxiliary data vector phi  are 
  21. %    obtained by [THM,YHAT,P,phi]=rarx(Z,NN,adm,adg,TH0,P0,phi0). 
  22. %    
  23. %    See also RARMAX, ROE, RBJ, RPEM and RPLR.
  24.  
  25. %    L. Ljung 10-1-89
  26. %    Copyright (c) 1989-90 by the MathWorks, Inc.
  27. %    All Rights Reserved.
  28.  
  29.  
  30. [nz,ns]=size(z);
  31. if ns==1,if length(nn)~=1,error('For a time series nn should be a scalar nn=na!'),end,end
  32. if 2*ns-1~=length(nn),error('Incorrect number of orders specified in nn. nn=[na nb nk]'),end
  33. na=nn(1);if ns>1,nb=nn(2:ns);nk=nn(ns+1:2*ns-1);else nk=1;nb=0;end
  34. nu=1;
  35. if any(nk<1),error('Sorry, this routine requires nk>0; Shift input sequence if necessary!'),end
  36. d=na+sum(nb);
  37. nbm=nb+nk-1;ncbm=na+cumsum([0 nbm]);
  38. ii=[1:na+sum(nbm)];
  39. i=[1:na];
  40. for ku=1:ns-1,i=[i ncbm(ku)+nk(ku):ncbm(ku+1)];end 
  41.  
  42. dm=na+sum(nbm);
  43.  
  44. if nargin<7, phi=zeros(dm,1);end
  45. if nargin<6, p0=10000*eye(d);end
  46. if nargin<5, th0=eps*ones(d,1);end
  47. if isempty(phi),phi=zeros(dm,1);end
  48. if isempty(p0),p0=10000*eye(d);end
  49. if isempty(th0),th0=eps*ones(d,1);end
  50. if length(th0)~=d, error('The length of th0 must equal the number of estimated parameters!'),end
  51. [th0nr,th0nc]=size(th0);if th0nr<th0nc, th0=th0';end
  52. p=p0;th=th0;
  53. if adm(1)=='f', R1=zeros(d,d);lam=adg;end
  54. if adm(1)=='k', [sR1,SR1]=size(adg);
  55.      if sR1~=d | SR1~=d,error('The R1 matrix should be a square matrix with dimension equal to number of parameters!'),end
  56.      R1=adg;lam=1;
  57. end
  58. if adm(2)=='g', grad=1;else grad=0;end
  59.  
  60. for kcou=1:nz
  61. yh=phi(i)'*th;
  62. epsi=z(kcou,1)-yh;
  63. if ~grad,K=p*phi(i)/(lam + phi(i)'*p*phi(i));
  64.          p=(p-K*phi(i)'*p)/lam+R1;
  65. else K=adg*phi(i);end
  66. if adm(1)=='n', K=K/(eps+phi(i)'*phi(i));end
  67. th=th+K*epsi;
  68.  
  69. epsilon=z(kcou,1)-th'*phi(i);
  70.  
  71. phi(ii+1)=phi(ii);
  72. if na>0,phi(1)=-z(kcou,1);end
  73. if any(ncbm>0),phi(ncbm(1:ns-1)+1)=z(kcou,2:ns)';end
  74.  
  75. thm(kcou,:)=th';yhat(kcou)=yh;
  76. end
  77. yhat=yhat';