home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / IDSIMSD.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  1.7 KB  |  58 lines
  1. function idsimsd(u,th,n,noise,ky)
  2. %IDSIMSD Illustrates the uncertainty in simulated model responses
  3. %
  4. %    idsimsd(U,TH)
  5. %
  6. %    U is a column vector (matrix) containing the input(s).
  7. %    TH is a model given in the THETA format (See HELP THETA).
  8. %    10 random models are created, consistent with the covariance informa-
  9. %    tion in TH, and the responses of each of these models to U are plotted
  10. %    in the same diagram.
  11. %    
  12. %    The number 10 can be changed to N by idsimsd(U,TH,N).
  13. %
  14. %    With idsimsd(U,TH,N,'noise',KY), additive noise (e) is added to the 
  15. %    simulation in accordance with the noise model of TH.
  16. %    KY denotes the output numbers to be plotted (default all)
  17.  
  18. %    L.Ljung 7-8-87
  19. %    Copyright (c) The MathWorks, Inc.
  20. %    All Rights Reserved
  21.  
  22. if isthss(th),ny=th(1,4);else ny=1;end
  23. if nargin<5,ky=[];end
  24. if nargin <4,noise=[];end
  25. if nargin<3,n=[];end,
  26. if isempty(ky),ky=1:ny;end
  27. if isempty(noise),noise='nonoise';end
  28. if isempty(n),n=10;end
  29.  
  30. nu=th(1,3);Tsamp=abs(th(1,2));
  31. [Nc,d]=getncap(th);
  32. [N,nz]=size(u);
  33. if nu~=nz, error('The input matrix U has not the correct number of columns!'),
  34. return,end
  35. [par,P,lam]=th2par(th);
  36. if size(P)==0; disp('No covariance information given in TH'),P=zeros(d,d);end
  37. P=chol(P);
  38. u1=u;
  39. for kk=ky
  40. if noise(3)=='i', u1=[u randn(N,ny)];end %corr 9007
  41. yh=idsim(u1,th);yh=yh(:,kk);
  42. ndu=length(yh);y1=max(yh);y2=min(yh);
  43. y12=y1-y2;y1=y1+0.2*y12;y2=y2-0.2*y12;
  44. plot([1:ndu]*Tsamp,yh),axis([Tsamp ndu*Tsamp y2 y1]);hold on;
  45. title(['Output number ',int2str(kk)])
  46. u1=u;
  47. for k=1:n-1
  48.     th1=th; th1(3,1:d)=par+randn(1,d)*P;
  49.     if noise(3)=='i', u1=[u randn(N,ny)];end
  50.     yh=idsim(u1,th1);
  51.     plot([1:ndu]*Tsamp,yh(:,kk))
  52. end
  53. hold off
  54. if kk<ky(length(ky)),pause,end
  55. end
  56. set(gcf,'NextPlot','replace');
  57.  
  58.