home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / SELSTRUC.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.7 KB  |  74 lines
  1. function [nn,Vmod]=selstruc(V,c)
  2. %SELSTRUC    Selects model structures according to various criteria
  3. %
  4. %    nn = selstruc(V,c)    or   [nn,Vm] = selstruc(V,c)
  5. %
  6. %    V: is a matrix containing information about different structures,
  7. %    typically obtained as the output of ARXSTRUC or IVSTRUC.
  8. %    
  9. %    c: selects the criterion: c='PLOT' gives plots of the loss function
  10. %    and (if applicable) the conditioning number from IVSTRUC as func-
  11. %    tions of the number of estimated parameters. The user then selects the
  12. %    number of parameters. 
  13. %    Automatic choices of structures are obtained by c='AIC', which gives
  14. %    Akaike's information theoretic criterion, while c='MDL' gives 
  15. %    Rissanen's minimum description length criterion. If c is given a 
  16. %    numeric value the structure is selected by minimization of
  17. %    (1 + c*d/N)*Vd, where d is the number of estimated parameters,
  18. %    Vd is the loss function of the corresponding model, and N is the
  19. %    number of data. If c is omitted, PLOT is chosen.
  20. %
  21. %    nn: is returned as the chosen structure. The format is compatible
  22. %    with the input format for ARX and IV4.
  23. %    Vm: the first row of Vm contains the logarithms of the modified
  24. %    criteria of fit. The remaining rows of Vm coincide with V.
  25.  
  26.  
  27. %    L. Ljung 4-12-87
  28. %    Copyright (c) 1987-92 by the MathWorks, Inc.
  29. %    All Rights Reserved.
  30.  
  31. if nargin<2,c='p';,end
  32. if c<0,c='p';end
  33. [nl1,nm1]=size(V);
  34. nu=floor((nl1-2)/2);
  35. Nc=V(1,nm1);
  36. if c(1)=='a' | c(1)=='A',alpha=2;end
  37. if c(1)=='m' | c(1)=='M', alpha=log(Nc);end
  38. if norm(c(1)=='mMaApPlL')==0,alpha=c;end
  39. if norm(c(1)=='PpLl')>0,clg,alpha=0;
  40.     if nu>0,sv=sum(V(2:2+nu,1:nm1-1));else sv=V(2,1:nm1-1);end %Cor 9007
  41.     vv=V(1,1:nm1-1);
  42.     if c(1)=='l' |c(1)=='L', vv=log(vv);end
  43.     kt=((2+2*nu)==nl1);
  44.     if kt==1,subplot(111),else subplot(121),end
  45.     axsv=[min(sv)-1 max(sv)+1];
  46.     plot(sv,vv,'*',axsv,[max(vv) min(vv)],'.')
  47.     xlabel('# of par`s'),
  48.     if c(1)=='l' |c(1)=='L',ylabel('log of loss fcn'),
  49.     else ylabel('loss fcn'),end
  50.     if kt==0, subplot(122),plot(sv,V(nl1,1:nm1-1),'o',axsv,[min(V(nl1,1:nm1-1)) max(V(nl1,1:nm1-1))],'.')
  51.         xlabel('# of par`s'),ylabel('log of cond`ing #'),end
  52.     title('RETURN TO COMMAND SCREEN TO SELECT # OF PARAMETERS TO BE ESTIMATED')
  53.     pause
  54.     if exist('testall')
  55.         pp=-1;
  56.     else    
  57.     pp=input('Enter # of parameters to be estimated ("enter" gives default): ');
  58.     end
  59.     if length(pp)>0 & pp>0,kk=find(sv==pp);
  60.     if length(kk)==0, disp('This number of par`s is not available'),return,end
  61.     [vm,skl]=min(V(1,kk));nn=V(2:2+2*nu,kk(skl))';
  62.     end
  63. end    
  64. for kj=1:nm1-1
  65. Vmod(1,kj)=V(1,kj)*(1+(alpha/Nc)*sum(V(2:nu+2,kj)));
  66. end
  67. Vmod(2:nl1,1:nm1-1)=V(2:nl1,1:nm1-1);
  68.  
  69. Vmod(1,1:nm1-1)=log(Vmod(1,1:nm1-1));
  70. if length(nn)==0,[vm,sel]=min(Vmod(1,1:nm1-1));
  71. nn=V(2:2+2*nu,sel)';end
  72. set(gcf,'NextPlot','replace');
  73.  
  74.