home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / SEGMENT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  4.0 KB  |  125 lines
  1. function [segm,V,thm,r2e]=segment(z,nn,r2,q,r1,M,th0,p0,lifelength,mu)
  2. %SEGMENT  Segments data and tracks abruptly changing systems 
  3. %
  4. %       [segm,V] = segment(z,nn,r2,q)
  5. %
  6. %       z : The output-input data vector z = [y u].
  7. %       nn: ARX or ARMAX model orders nn = [na nb nk] or nn = [na nb nc nk].
  8. %        See HELP ARX or ARMAX. The algorithm handles multi-input systems.
  9. %       r2: The equation error variance(Default:estimated, but better to guess)
  10. %    q : The probability that the system jumps at each sample.(Default 0.01)
  11. %       segm: The parameters of the segmented data. Row k is for sample # k
  12. %          The parameters are given in "alphabetical order". 
  13. %    V: The loss function corresponding to segm
  14. %    The time-varying estimates th, and the estimated values of r2 are given
  15. %    by   [segm,V,th,r2] = segment(z,nn)
  16. %    
  17. %    Additional design variables are reached by
  18. %    [segm,V,th,r2]=segment(z,nn,r2,q,R1,M,th0,P0,ll,mu)
  19. %    R1: covariance matrix of the parameter jumps. (Default Identity matrix)
  20. %    M: Number of parallel models to be used (Default 5)
  21. %    th0:Initial parameter estimates (row vector)
  22. %    P0:Initial covariance matrix (Def 10*I). ll:Guaranteed lifelength of
  23. %    each model(Def 1). mu: Forgetting factor in r2-estimation (Def 0.97).
  24. %    Reference: P. Andersson Int J Control, Nov 1985.
  25.  
  26. %    L. Ljung 10-1-89
  27. %    Copyright (c) 1989 by the MathWorks, Inc.
  28. %    All Rights Reserved.
  29.  
  30. [Ncap,nz]=size(z);
  31. if Ncap<nz,disp('Warning: Have you entered the data as column vectors?'),end
  32. nu=nz-1; nl=length(nn);
  33. if nu==0,na=nn(1);nb=0;nc=0;nk=1;if nl>1,nc=nn(2);end
  34. if nl>2,error('For a time series, nn = na or nn = [na nc]!'),end
  35. else
  36.   na=nn(1);nb=nn(2:nu+1);
  37.   if nl==2*nu+1,nk=nn(nu+2:2*nu+1);nc=0;
  38.   else if nl==2*nu+2,nc=nn(nu+2);nk=nn(nu+3:2*nu+2);
  39.   else error('Incorrect number of orders specified: nn=[na nb nk] or nn=[na nb nc nk]'),end,end
  40. end
  41. d=na+sum(nb)+nc; % Number of parameters
  42. if nargin<10, mu=[];end
  43. if nargin<9, lifelength=[];end
  44. if nargin<8 p0=[];end
  45. if nargin<7,th0=[];end
  46. if nargin<6 M=[];end
  47. if nargin<5 r1=[];end
  48. if nargin<4 q=[];end
  49. if nargin<3 r2=[];end
  50. if isempty(q),q=0.01;end
  51. if isempty(r1),r1=eye(d);end
  52. if isempty(M),M=5;end
  53. if isempty(th0),th0=zeros(1,d);end
  54. if isempty(p0),p0=10*eye(d);end
  55. if isempty(lifelength),lifelength=1;end
  56. if isempty(mu),mu=0.97;end
  57. if isempty(r2),r2dum=1;R2v=ones(1,M);r2=r2dum;estr2=1;else estr2=0; end;
  58. alfa=1/M*ones(1,M);%[1,zeros(1,M-1)];
  59.  
  60. thm=zeros(Ncap,d);
  61. Phicap=zeros(Ncap,d);
  62. th=th0'*ones(1,M);age=zeros(1,M);hist=zeros(Ncap,M);r2e=zeros(Ncap,1);
  63. P=[]; for i=1:M,P=[P p0/r2];end;
  64.  
  65. if nc>0,nu=nu+1;nk(nu)=1;nb(nu)=nc;z(1:nc,nu+1)=zeros(nc,1);end
  66. for i=max([na+1,max(nb+nk),nc+1]):Ncap
  67.     phi=[-z(i-1:-1:i-na,1)];
  68.     for ku=1:nu,
  69.     phi=[phi;z(i-nk(ku):-1:i-nk(ku)-nb(ku)+1,1+ku)];
  70.     end
  71.     y=z(i,1);
  72.     for j=1:M,
  73.         Pj=P(:,d*(j-1)+1:d*(j-1)+d);
  74.         den(j)=(r2+phi'*Pj*phi);
  75.         
  76.         epsi(j)=y-th(:,j)'*phi;
  77.  
  78.         alfabar(j)=alfa(j)/sqrt(den(j))*exp(-epsi(j)^2/(2*den(j)));
  79.         th(:,j)=th(:,j)+1/den(j)*P(:,d*(j-1)+1:d*(j-1)+d)*phi*epsi(j);
  80.     P(:,d*(j-1)+1:d*(j-1)+d)=Pj-(Pj*phi*phi'*Pj)/den(j);
  81.     end;
  82.  
  83.     aind=find((age(1:M)>lifelength));if length(aind)>0,
  84.     [dummy,jmin]=min(alfabar(aind));jmin=aind(jmin);
  85.     [dummy,jmax]=max(alfabar);
  86.  
  87.     P(:,d*(jmin-1)+1:d*(jmin-1)+d)=P(:,d*(jmax-1)+1:d*(jmax-1)+d)+r1;
  88.     th(:,jmin)=th(:,jmax);
  89.     alfabar(jmin)=q*alfabar(jmax);age(jmin)=0;
  90.     hist(:,jmin)=hist(:,jmax);hist(i,jmin)=1;
  91.     else jmax=1;end
  92.     age=age+ones(1,M);
  93.     alfa=1/sum(alfabar)*alfabar;
  94.  
  95.     theta=th*alfa';
  96.     epsilon=y-theta'*phi;if nc>0,z(i,nu+1)=epsilon;end
  97. %
  98. %    estimate R2
  99. %
  100. if estr2
  101.     r2dum=r2dum+(1-mu)*(epsilon^2-r2dum);
  102.     agedum=max(age,2)-1;
  103.     R2v= R2v +max(ones(1,M)./agedum,(1-mu)).*(epsi.^2./den-R2v); 
  104.     Rr=R2v(find(age>d));
  105.     if length(Rr)>0,r2=min(Rr);else r2=r2dum;end
  106.  
  107. end
  108.  
  109. %
  110.  
  111. %
  112.     thm(i,:)=theta';r2e(i)=r2;Phicap(i,:)=phi';
  113.  
  114. %
  115. end
  116. hh=hist(:,jmax);
  117. kk=find(hh==1);
  118. kk=[1 kk' Ncap];
  119. nn=length(kk);
  120. for k=2:nn
  121. segm(kk(k-1):kk(k),:)=ones(kk(k)-kk(k-1)+1,1)*thm(kk(k),:);
  122. end
  123. e=z(:,1)-sum(segm'.*Phicap')';
  124. V=e'*e/Ncap;
  125.