home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / RPEM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  5.0 KB  |  141 lines
  1. function [thm,yhat,p,phi,psi] = rpem(z,nn,adm,adg,th0,p0,phi,psi)
  2. %RPEM    Computes estimates recursively for a general model using the
  3. %    Recursive Prediction Error Method
  4. %
  5. %    [THM,YHAT] = rpem(Z,NN,adm,adg)
  6. %
  7. %    Z: The output-input data z=[y u] (Multiple input is allowed)
  8. %    NN : NN=[na nb nc nd nf nk], The orders and delay of a general
  9. %         input-output model (see HELP PEM)
  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] = rpem(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]=rpem(Z,NN,adm,adg,TH0,P0,phi0,psi0). 
  23. %    
  24. %    See also RARMAX, RARX, ROE, RBJ 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);nu=ns-1;
  31. if ns==1,error('Use RARMAX instead for a time series!'),end
  32. nll=length(nn);
  33. if length(nn)~=3+3*nu,error('An incorrect number of orders is specified. nn should be [na nb nc nd nf nk] with nb, nf and nk of length = # of inputs!'),end
  34. na=nn(1);nb=nn(2:ns);nc=nn(ns+1);nd=nn(ns+2);
  35. if any(nb==0),error('If (some) nb is zero, please remove that input first!'),end
  36. nf=nn(ns+3:ns+2+nu);nk=nn(ns+3+nu:ns+2+2*nu);
  37. if any(nk<1),error('Sorry, this routine requires nk>0; Shift input sequence if necessary!'),end
  38. d=sum(nn(1:ns+2+nu));
  39. ng=nc+nf;
  40. nam=(na>0)*max([na,nd,nc]);
  41. nbm=(nb>0).*max(max(nb,nd)+nk-1,ng);
  42. ncm=nc;
  43. ndm=(nd>0)*max(nd,nc);
  44. nfm=(nf>0).*max(nd,ng);
  45.  
  46. %First indices for each par:
  47. ina=(na>0)*1;
  48. inbt=(nam+cumsum([0 nbm])+1);
  49. inb=(nb>0).*inbt(1:nu);
  50. inct=inbt(ns);
  51. inc=(nc>0)*inct;
  52. indt=(inct+ncm);
  53. ind=(nd>0)*indt;
  54. infft=(indt+ndm+cumsum([0 nfm]));
  55. inff=(nf>0).*infft(1:nu);
  56. lastind=infft(nu)+nfm(nu)-1;
  57. initind=[ina inb inc ind inff];
  58.  
  59.  %Indices in theta-vector
  60. tia=1:na;tib=na+1:na+sum(nb);
  61. tic=na+sum(nb)+1:na+sum(nb)+nc;
  62. nnid=na+sum(nb)+nc+nd;tid=na+sum(nb)+nc+1:nnid;
  63. tif=nnid+1:nnid+sum(nf);
  64.  
  65. ia=1:na;
  66. ic=inc:inc+nc-1;
  67. id=ind:ind+nd-1;
  68. ibs=[];ifs=[];for ku=1:nu,ibs=[ibs inb(ku)-1+nk(ku):inb(ku)+nk(ku)+nb(ku)-2];
  69. ifs=[ifs inff(ku):inff(ku)+nf(ku)-1];end
  70. i=[ia ibs ic id ifs];
  71. dm=lastind;
  72.  
  73. if nargin<8, psi=zeros(dm,1);end
  74. if nargin<7, phi=zeros(dm,1);end
  75. if nargin<6, p0=10000*eye(d);end
  76. if nargin<5, th0=eps*ones(d,1);end
  77. if isempty(psi),psi=zeros(dm,1);end
  78. if isempty(phi),phi=zeros(dm,1);end
  79. if isempty(p0),p0=10000*eye(d);end
  80. if isempty(th0),th0=eps*ones(d,1);end
  81. if length(th0)~=d,error('The length of th0 must equal the number of estimated parameters!'),end
  82. [th0nr,th0nc]=size(th0);if th0nr<th0nc,th0=th0';end
  83. p=p0;th=th0;
  84. if adm(1)=='f', R1=zeros(d,d);lam=adg;end
  85. if adm(1)=='k', [sR1,SR1]=size(adg);
  86.      if sR1~=d | SR1~=d,error('The R1 matrix should be a square matrix with dimension equal to number of parameters!'),end
  87.      R1=adg;lam=1;
  88. end
  89. if adm(2)=='g', grad=1;else grad=0;end
  90.  
  91. for kcou=1:nz
  92. yh=phi(i)'*th;
  93. epsi=z(kcou,1)-yh;
  94. if ~grad,K=p*psi(i)/(lam + psi(i)'*p*psi(i));
  95.          p=(p-K*psi(i)'*p)/lam+R1;
  96. else K=adg*psi(i);end
  97. if adm(1)=='n', K=K/(eps+psi(i)'*psi(i));end
  98. th=th+K*epsi;
  99. d=[1;th(tid)];
  100. if nc>0,c=fstab([1;th(tic)])';else c=1;end,th(tic)=c(2:nc+1);
  101. if ~isempty(tif),tif3=tif(1);else tif3=0;end
  102. if ~isempty(tib),tib3=tib(1);else tib3=0;end
  103. for ku=1:nu,
  104.     if nf(ku)>0,f=fstab([1;th(tif3:tif3+nf(ku)-1)])';else f=1;end,g=conv(f,c);
  105.     th(tif3:tif3+nf(ku)-1)=f(2:nf(ku)+1);
  106.     if nb(ku)+nf(ku)>0,w(ku)=th([tib3:tib3+nb(ku)-1 tif3:tif3+nf(ku)-1])'*...
  107.        phi([inb(ku)+nk(ku)-1:inb(ku)+nk(ku)+nb(ku)-2 inff(ku):inff(ku)+nf(ku)-1]);
  108.     else w(ku)=0;
  109.     end 
  110.     tif3=tif3+nf(ku);tib3=tib3+nb(ku);
  111.     if nb(ku)>0,
  112.        util(ku)=d'*[z(kcou,ku+1);phi(inb(ku)+nk(ku)-1:inb(ku)+nk(ku)+nd-2)]-...
  113.        g'*[0;psi(inb(ku):inb(ku)+ng(ku)-1)];
  114.     end
  115.     if nf(ku)>0
  116.        wtil(ku)=d'*[w(ku);-phi(inff(ku):inff(ku)+nd-1)]+g'*[0;psi(inff(ku):inff(ku)+ng(ku)-1)];
  117.     end
  118. end
  119. v=[z(kcou,1);-phi(ia)]'*[1;th(ia)]-sum(w);
  120. epsilon=v-th([tic tid])'*phi([ic id]);
  121.  
  122. if na>0,ytil=d'*[z(kcou,1);-phi(1:nd)]+c'*[0;psi(1:nc)];end
  123.  
  124.  
  125. if nc>0
  126. epstil=c'*[epsilon;-psi(ic)];end
  127. if nd>0
  128. vtil=c'*[v;psi(ind:ind+nc-1)];
  129. end
  130. phi(2:dm)=phi(1:dm-1);psi(2:dm)=psi(1:dm-1);
  131. if na>0,phi(1)=-z(kcou,1);psi(1)=-ytil;end
  132. for ku=1:nu
  133.     if nb(ku)>0,phi(inb(ku))=z(kcou,1+ku);psi(inb(ku))=util(ku);end 
  134.     if nf(ku)>0,phi(inff(ku))=-w(ku);psi(inff(ku))=-wtil(ku);end
  135. end
  136. if nc>0,phi(inc)=epsilon;psi(inc)=epstil;end
  137. if nd>0,phi(ind)=-v;psi(ind)=-vtil;end
  138.  
  139. thm(kcou,:)=th';yhat(kcou)=yh;
  140. end
  141. yhat=yhat';