home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / MINCRIT.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.0 KB  |  69 lines
  1. function THcap=mincrit(z,th,index,maxiter,tol,lim,maxsize,T)
  2. %MINCRIT   Minimizes prediction error criteria
  3. %
  4. %    TH = mincrit(Z,TH,MAXITER,TOL,LIM,MAXSIZE,T)
  5. %
  6. %    This routine minimized robustified quadratic prediction error
  7. %    criteria. It is intended primarily as a subroutine to PEM.
  8. %    See PEM for an explanation of arguments.
  9.  
  10. %    L. Ljung 10-1-86
  11. %    Copyright (c) 1986-92 by the MathWorks, Inc.
  12. %    All Rights Reserved.
  13.  
  14. nu=th(1,3);
  15. NN=th(1,4:6+2*nu);n=sum(NN); [Ncap,nz]=size(z);
  16. if isempty(index),index=1:n;end
  17. g=ones(n,1);, l=0;
  18. e=pe(z,th);
  19.  
  20.  
  21. %** Determine limit for robustification of criterion ***
  22.  
  23. if lim~=0, lim=median(abs(e-median(e)))*lim/0.7;end
  24. if lim==0, lim=1000*th(1,1);end
  25.  
  26. [me,ne]=size(e);
  27. el=min(e,lim*ones(me,ne)); el=max(el,-lim*ones(me,ne)); th(1,1)=el'*e/length(e);
  28. clear e, l=0; st=0;
  29.  
  30. %  ** The minimization loop **
  31.  
  32. while [norm(g)>tol l<maxiter st~=1]
  33.    l=l+1;
  34.  
  35. %   * Compute the Gauss-Newton direction *
  36.  
  37.    [g,R,grad]=gn(z,th,lim,maxsize,index);
  38.  
  39. %    * Search along the g-direction for a lower value of V *
  40.  
  41.     [th1,st]=search(z,th,g,lim);
  42.     if st==1,
  43.     [th1,st]=search(z,th,grad/trace(R)*length(R),lim);end
  44. %  * Display the current values along with the previous ones *
  45.  
  46.      t=zeros(n,3); t(:,1)=th1(3,1:n)';t(:,2)=th(3,1:n)';t(:,3)=g;
  47.      home
  48.      disp(['    ITERATION # ' int2str(l)])
  49.      disp(['Current loss: ' num2str(th1(1,1)) '   Previous loss: ' num2str(th(1,1))])
  50.       disp(['Current th  prev. th   gn-dir. '])
  51.       disp(t)
  52.       disp(['Norm of gn-vector: ' num2str(norm(g))])
  53.       if st==1
  54.           disp('No improvement of the criterion possible along the search direction')
  55.           disp('Iterations therefore terminated')
  56.       end
  57.       th=th1;
  58. end
  59.  
  60. % *** Build up the TH-matrix ***
  61.  
  62. THcap=th;THcap(1,1)=th(2,1);THcap(1,2)=T;
  63. THcap(2,1)=(1+n/Ncap)/(1-n/Ncap)*THcap(1,1);
  64. THcap(4:n+3,1:n)=zeros(n);
  65. THcap(index+3,index)=THcap(1,1)*inv(R);
  66. ti=fix(clock); ti(1)=ti(1)/100;
  67. THcap(2,2:6)=ti(1:5); THcap(2,7)=5;
  68.  
  69.