home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / ARXSTRUC.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  3.3 KB  |  100 lines
  1. function V=arxstruc(z,zv,nn,maxsize);
  2. %ARXSTRUC    Computes loss functions for families of ARX-models
  3. %
  4. %    V = arxstruc(Z,ZV,NN)
  5. %
  6. %    V: The first row of V is returned as the loss functions of the 
  7. %    structures defined by the rows of NN. The remaining rows are 
  8. %    returned as NN'. The last column of V contains the number of 
  9. %    data points in Z.
  10. %
  11. %    Z : the output-input data Z=[y u], with y and u as column vectors.
  12. %    Only Single-Output data are handled by this routine!
  13. %    For a time series Z=y only.
  14. %    ZV: the input-output data on which the validation is performed.
  15. %    (could equal Z). The number of rows in ZV need not equal those in Z.
  16. %
  17. %    NN: Each row of NN is of the format [na nb nk], the orders and delays
  18. %    of the corresponding ARX model.(See HELP ARX and HELP STRUC.)
  19. %    For a time series, each row of NN equals na only.
  20. %
  21. %    See HELP SELSTRUC for analysis of V.
  22. %    Some parameters associated with the algorithm are accessed by
  23. %    V = arxstruc(Z,ZV,NN,maxsize)
  24. %    See HELP AUXVAR for an explanation of maxsize.
  25.  
  26. %    
  27. %    L. Ljung 4-12-87
  28. %    Copyright (c) 1987-90 by the MathWorks, Inc.
  29. %    All Rights Reserved.
  30.  
  31. [Ncap,nz]=size(z); nu=nz-1;[nm,nl]=size(nn);[Ncapv,nzv]=size(zv);
  32. if nu>1 & nm==1,error('For just one model, use ARX instead!'),end
  33. if nz>Ncap,error(' Data should be organized in column vectors!'),return,end
  34. if nl~=1+2*(nz-1),disp(' Incorrect number of orders specified!'),disp('For an AR-model nn=na'),disp('For an ARX-model, nn=[na nb nk]'),
  35. error('see above'),return,end
  36. if nzv~=nz,error(' ZV-set should have the same number of inputs as Z!'),return,end
  37. if nz>1, na=nn(:,1);nb=nn(:,2:1+nu);nk=nn(:,2+nu:1+2*nu);
  38. else na=nn(:,1); nb=zeros(nm,1);nk=zeros(nm,1);end
  39. nma=max(na);nbkm=max(nb+nk)-ones(1,nu);nkm=min(nk);
  40. n=nma+sum((nbkm-nkm))+nu;
  41. if nn==0,return,end
  42. vz=1;if size(z)==size(zv),if norm(z-zv)<eps, vz=0;end,end
  43.  
  44. % *** Set up default values **
  45. maxsdef=idmsize(Ncap,n);
  46.  
  47. if nargin<4, maxsize=maxsdef;end
  48. if isempty(maxsize),maxsize=maxsdef;end
  49. if maxsize<0, maxsize=maxsdef; end
  50.  
  51. % *** construct regression matrix ***
  52.  
  53.  nmax=max(max([na+ones(nm,1) nb+nk]'))-1;
  54. M=floor(maxsize/n);
  55. R=zeros(n);F=zeros(n,1);if vz,Rv=zeros(n);Fv=zeros(n,1);end
  56. for k=nmax:M:max(Ncap,Ncapv)
  57.     if min(Ncap,k+M)<k+1,ntz=0;else ntz=1;end
  58.     if min(Ncapv,k+M)<k+1,ntzv=0;else ntzv=1;end
  59.       if ntz,jj=(k+1:min(Ncap,k+M));phi=zeros(length(jj),n);end
  60.       if vz & ntzv,jjv=(k+1:min(Ncapv,k+M));phiv=zeros(length(jjv),n);end
  61.       for kl=1:nma, 
  62.       if ntz,phi(:,kl)=-z(jj-kl,1);end
  63.       if vz & ntzv ,phiv(:,kl)=-zv(jjv-kl,1);end,end
  64.       ss=nma;
  65.       for ku=1:nu
  66.            for kl=nkm(ku):nbkm(ku),
  67.            if ntz,phi(:,ss+kl+1-nkm(ku))=z(jj-kl,ku+1);end
  68.            if vz & ntzv,phiv(:,ss+kl+1-nkm(ku))=zv(jjv-kl,ku+1);end
  69.        end    
  70.            ss=ss+nbkm(ku)-nkm(ku)+1;
  71.       end
  72.       if ntz,R=R+phi'*phi; F=F+phi'*z(jj,1);end
  73.       if vz & ntzv,Rv=Rv+phiv'*phiv; Fv=Fv+phiv'*zv(jjv,1);end
  74. end
  75. %
  76. % *** compute loss functions ***
  77. %
  78. v1=zv(nmax+1:Ncapv,1)'*zv(nmax+1:Ncapv,1);
  79. for j=1:nm
  80.     s=[1:na(j)];rs=nma;
  81.     for ku=1:nu
  82.         
  83.         s=[s,rs+nk(j,ku)-nkm(ku)+1:rs+nb(j,ku)+nk(j,ku)-nkm(ku)];
  84.         rs=rs+nbkm(ku)-nkm(ku)+1;    
  85.     end
  86.     
  87.     RR=R(s,s);
  88.     FF=F(s);
  89.     if vz,RRv=Rv(s,s);FFv=Fv(s);end
  90.     TH=(RR\FF);
  91.     if vz,
  92.         V(1,j)=(v1-2*FFv'*TH + TH'*RRv*TH)/(Ncapv-nmax);
  93.     else    
  94.         V(1,j)=(v1-FF'*TH)/(Ncap-nmax);
  95.     end
  96.     V(2:nl+1,j)=nn(j,:)';
  97.  
  98. end;
  99. V(1,nm+1)=Ncap-nmax;
  100.