home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / PEM.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.9 KB  |  93 lines
  1. function th=pem(z,nn,index,maxiter,tol,lim,maxsize,Tsamp)
  2. %PEM    Computes the prediction error estimate of a general linear model
  3. %    TH = pem(Z,THSTRUC)
  4. %
  5. %    TH : returned as the estimated parameters of the model
  6. %    along with estimated covariances and structure information. 
  7. %    For the exact format of TH see HELP THETA.
  8. %
  9. %    Z : the output input data [y u], with y and u as column vectors
  10. %    For multi-variable systems, Z=[y1 y2 ... yp u1 u2 ... un]
  11. %
  12. %     THSTRUC: A matrix that defines the model structure. Typically created
  13. %    by POLY2TH, MS2TH or MF2TH or by some estimation routine. The minimiza-
  14. %    tion is then initialized at the parameters contained in THSTRUC. 
  15. %    A general, black-box multi-input structure 
  16. %    A(q) y(t) = [B(q)/F(q)] u(t-nk) + [C(q)/D(q)] e(t)
  17. %    is obtained for    THSTRUC = [na nb nc nd nf nk], indicating the orders
  18. %    in the above model. 
  19. %    By TH=pem(Z,THSTRUC,INDEX) only parameters corresponding to the indices
  20. %    in the row vector INDEX are estimated.
  21. %    Some parameters associated with the algorithm are accessed by
  22. %    TH = pem(Z,THSTRUC,index,maxiter,tol,lim,maxsize,T)
  23. %    See HELP AUXVAR for an explanation of these and their default values.
  24.  
  25. %    L. Ljung 10-1-86, 10-2-90
  26. %    Copyright (c) 1986-92 by the MathWorks, Inc.
  27. %    All Rights Reserved.
  28.  
  29. %  *** Set up default values ***
  30.  
  31. [Ncap,nz]=size(z);
  32. maxsdef=idmsize(Ncap);
  33.  
  34. if nargin<8, Tsamp=[];end
  35. if nargin<7, maxsize=[];end
  36. if nargin<6, lim=[];end
  37. if nargin<5, tol=[];end
  38. if nargin<4, maxiter=[];end
  39. if nargin<3, index=[];end
  40. [nr,cc]=size(nn);
  41. if isempty(Tsamp),if nr>1,Tsamp=nn(1,2);else Tsamp=1;end,end
  42. if isempty(maxsize),maxsize=maxsdef;end,if maxsize<0,maxsize=maxsdef;end
  43. if isempty(lim),lim=1.6;end,if lim<0,lim=1.6;end
  44. if isempty(tol),tol=0.01;end,if tol<0,tol=0.01;end
  45. if isempty(maxiter),maxiter=10;end,if maxiter<0,maxiter=10;end
  46. if nr>1,if isthss(nn)
  47.     eval('th=pemss(z,nn,index,maxiter,tol,lim,maxsize,Tsamp);')
  48. return,end,end 
  49.  
  50. if Tsamp<0, Tsamp=1;end
  51.  
  52. % *** Do some consistency tests ***
  53.  
  54. nu=nz-1;
  55. if nz>Ncap, error('The data should be organized in column vectors!')
  56. return,end
  57. nb=nn(2:nz); if nz==1 
  58. error('Use ARMAX instead for a time series!')
  59. end
  60. if nr==1 & cc~=3+3*nu,
  61.     disp('Incorrect number of orders specified:')
  62.     disp('You should have nn=[na nb nc nd nf nk],')
  63.     disp('where nb, nf and nk are row vectors of length equal to')
  64.     disp(' the number of inputs'),error('see above')
  65. return,end
  66.  
  67.  
  68. %*** if nn has more than 1 row, then it is a theta matrix
  69. %    and we jump directly to mincrit ***
  70. %
  71. if nr>1,th=nn;end
  72.  
  73.  
  74. if nr==1
  75.    th=inival(z,nn,maxsize,Tsamp);
  76. end
  77.  
  78. % *** Display the initial estimate ***
  79.  
  80. clc, disp([' INITIAL ESTIMATE'])
  81. disp(['Current loss:' num2str(th(1,1))])
  82. disp(['th-vector'])
  83. nu=th(1,3); n=sum(th(1,4:6+2*nu));
  84. disp(th(3,1:n)')
  85.  
  86. if maxiter==0,return,end
  87.  
  88. % *** Minimize the prediction error criterion ***
  89.  
  90. th=mincrit(z,th,index,maxiter,tol,lim,maxsize,Tsamp);
  91.  
  92.  
  93.