home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / GN.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.9 KB  |  63 lines
  1. function [g,R,grad]=gn(z,th,lim,maxsize,index)
  2. %GN    Computes the Gauss-Newton direction for PE-criteria
  3. %
  4. %    [G,R] = gn(Z,TH,LIM,MAXSIZE,INDEX)
  5. %
  6. %    G : The Gauss-Newton direction
  7. %    R : The Hessian of the criterion (using the GN-approximation)
  8. %
  9. %    The routine is intended primarily as a subroutine to MINCRIT.
  10. %    See PEM for an explanation of arguments.
  11.  
  12. %    L. Ljung 10-1-86
  13. %    Copyright (c) 1986 by the MathWorks, Inc.
  14. %    All Rights Reserved.
  15.  
  16. % *** Set up the model orders ***
  17.  
  18. nu=th(1,3); na=th(1,4); nb=th(1,5:4+nu); nc=th(1,5+nu);
  19. nd=th(1,6+nu); nf=th(1,7+nu:6+2*nu); nk=th(1,7+2*nu:6+3*nu);
  20. [Ncap,nz]=size(z); n=na+sum(nb)+nc+nd+sum(nf);
  21. nmax=max([na nb+nk-ones(1,nu) nc nd nf]);
  22.  
  23. % *** Prepare for gradients calculations ***
  24.  
  25. [e,v,w,a,b,c,d,f]=pe(z,th);
  26. [me,ne]=size(e);
  27. ll=lim*ones(me,ne); el=min(e,ll); el=max(el,-ll); clear ll;
  28. if sum(nf)==0 , clear w, end
  29. if na>0, yf=filter(-d,c,z(:,1));end
  30. if nc>0, ef=filter(1,c,e); end
  31. if nd>0, vf=filter([-1],c,v); end
  32. for k=1:nu
  33.       gg=conv(c,f(k,:));
  34.       uf(:,k)=filter(d,gg,z(:,k+1));
  35.       if nf(k)>0, wf(:,k)=filter(-d,gg,w(:,k));end
  36. end
  37.  
  38. % *** Compute the gradient PSI. If N>M do it in portions ***
  39.  
  40. M=floor(maxsize/n);n1=length(index);
  41. R=zeros(n1);Fcap=zeros(n1,1);
  42. for k=nmax:M:Ncap-1
  43.       jj=(k+1:min(Ncap,k+M));
  44.       psi=zeros(length(jj),n);
  45.       for kl=1:na, psi(:,kl)=yf(jj-kl);,end
  46.       ss=na;ss1=na+sum(nb)+nc+nd;
  47.       for ku=1:nu
  48.              for kl=1:nb(ku), psi(:,ss+kl)=uf(jj-kl-nk(ku)+1,ku);end
  49.              for kl=1:nf(ku), psi(:,ss1+kl)=wf(jj-kl,ku);end
  50.              ss=ss+nb(ku);ss1=ss1+nf(ku);
  51.       end
  52.       for kl=1:nc, psi(:,ss+kl)=ef(jj-kl);end 
  53.       ss=ss+nc;
  54.       for kl=1:nd, psi(:,ss+kl)=vf(jj-kl);,end
  55.       psi=psi(:,index);
  56.       R=R+psi'*psi;  Fcap=Fcap+ psi'*el(jj);
  57. end
  58.  
  59. % *** Compute the Gauss-Newton direction g ***
  60.  
  61. if Ncap>M, g=R\Fcap; else g=psi\el(jj);end
  62. g1=zeros(n,1);g1(index)=g;g=g1;grad=zeros(n,1);grad(index)=Fcap;
  63.