Documentation of mincrit

Global Index (short | long) | Local contents | Local Index (short | long)

Function Synopsis

[THcap,it_inf]=mincrit(z,th,index,maxiter,tol,lim,maxsize,T,display,gui)

Help text

MINCRIT   Minimizes prediction error criteria

   TH = mincrit(Z,TH,MAXITER,TOL,LIM,MAXSIZE,T)

   This routine minimized robustified quadratic prediction error
   criteria. It is intended primarily as a subroutine to PEM.
   See PEM for an explanation of arguments.

Cross-Reference Information

Listing of function mincrit

function [THcap,it_inf]=mincrit(z,th,index,maxiter,tol,lim,maxsize,T,display,gui)

%   L. Ljung 10-1-86,9-25-93
%   Copyright (c) 1986-98 by The MathWorks, Inc.
%   $Revision: 2.4 $  $Date: 1997/12/03 23:16:14 $

global XIDiter
stopp=0;

nu=th(1,3);
NN=th(1,4:6+2*nu);n=sum(NN); [Ncap,nz]=size(z);
if isempty(index),index=1:n;end
if max(index)>n
  msg=str2mat('You have specified a parameter index number that is larger',...
                 'than the number of adjustable parameters in the model');
  if gui,errordlg(msg);iduipoin(2);end
  disp(msg),error(' ')
end
g=ones(n,1);, l=0;
e=pe(z,th);


%** Determine limit for robustification of criterion ***

if lim~=0, lim=median(abs(e-median(e)))*lim/0.7;end
if lim==0
    el=e;
else
    [ne,me]=size(e);
    ll=ones(ne,me)*lim;la=abs(e)+eps*ll;el=e.*(min(la,ll)./la);
    clear ll,clear la
end

th(1,1)=real(el'*e/length(e));
clear e, l=0; st=0;
if gui,drawnow,display=get(XIDiter(2),'value');end

%  ** The minimization loop **
testnorm=tol+1;
while [testnorm>tol l<maxiter st~=1]
   l=l+1;

%   * Compute the Gauss-Newton direction *

   [g,R,grad]=gn(z,th,lim,maxsize,index);
   testnorm=abs(grad'*g/th(1,1)/Ncap*100);

%    * Search along the g-direction for a lower value of V *

    [th1,st]=search(z,th,g,lim,display);
    if st==1,
    [th1,st]=search(z,th,grad/trace(R)*length(R),lim,display);end
%  * Display the current values along with the previous ones *

     t=zeros(n,3); t(:,1)=th1(3,1:n).';t(:,2)=th(3,1:n).';t(:,3)=g;
if isempty(t)|any(any(isinf(t)))
   disp('WARNING: Your input signal is apparently not persistently exciting.')
   disp('         Use other input or lower model orders.')
   disp('         No model returned.')
   th=[];return
end

if gui,drawnow,display=get(XIDiter(2),'value');end
if display
     %home
     disp(['    ITERATION # ' int2str(l)])
     disp(['Current loss: ' num2str(th1(1,1)) '   Previous loss: ' num2str(th(1,1))])
      disp(['Current th  prev. th   gn-dir. '])
      disp(t)
      disp(['Norm of gn-vector: ' num2str(norm(g))])
      disp(['Expected improvement: ',num2str(testnorm),' %'])
      disp(['Achieved improvement: ',num2str((th(1,1)-th1(1,1))*100/th(1,1)),' %'])

      if st==1
          disp('No improvement of the criterion possible along the search direction')
          disp('Iterations therefore terminated')
      end
end %if display
if gui
  set(XIDiter(4),'string',int2str(l));
  set(XIDiter(5),'string',num2str(th1(1,1)));
  set(XIDiter(9),'string',['(',num2str(th(1,1)),')']);
  set(XIDiter(6),'string',num2str(testnorm));
  drawnow
  stopp=get(XIDiter(7),'userdata');
end
      impr=th1(1,1)-th(1,1);th=th1;
      if stopp,break,end
end
it_inf=[l,impr,testnorm];
% *** Build up the TH-matrix ***

THcap=th;THcap(1,1)=th(2,1);THcap(1,2)=T;
THcap(2,1)=(1+n/Ncap)/(1-n/Ncap)*THcap(1,1);
THcap(4:n+3,1:n)=zeros(n);
THcap(index+3,index)=THcap(1,1)*inv(R);
ti=fix(clock); ti(1)=ti(1)/100;
THcap(2,2:6)=ti(1:5); THcap(2,7)=5;