home *** CD-ROM | disk | FTP | other *** search
/ Liren Large Software Subsidy 10 / 10.iso / l / l455 / 9.ddi / IDENT.DI$ / COMPARE.M < prev    next >
Encoding:
Text File  |  1993-03-11  |  2.0 KB  |  55 lines
  1. function [yh,fit] = compare(z,th,m,nr,za)
  2. % COMPARE  Compares the simulated/predicted output with the measured output
  3. %
  4. %    YH = compare(Z,TH,M)
  5. %
  6. %    Z : The output - input data for which the comparison is made (the
  7. %        validation data set.
  8. %    TH: The model in the THETA format (see help theta)
  9. %    M : The prediction horizon. Old outputs up to time t-M are used to
  10. %        predict the output at time t. All relevant inputs are used.
  11. %        M = inf gives a pure simulation of the system.(Default M=inf).
  12. %    YH: The resulting simulated/predicted output.
  13. %    COMPARE also plots YH together with the measured output in Z, and
  14. %        displays the mean square fit between these two signals.
  15. %        (solid/red is YH. dashed/green is measured output)
  16. %
  17. %    [YH,FIT] = compare(Z,TH,M,SAMPNR,LEVELADJUST)  
  18. %    gives access to some options:
  19. %    FIT: The mean square fit.
  20. %    SAMPNR: The sample numbers from Z to be plotted and used for the
  21. %       computation of FIT. (Default: SAMPNR = all rows of Z)
  22. %    LEVELADJUST: 'yes' adjusts the first values of YH and Z(:,1) to
  23. %      zero before plot and computation of FIT. 'no' is default.
  24.  
  25. %    L. Ljung 10-1-89,11-02-91
  26. %    Copyright (c) 1989-92 by the MathWorks, Inc.
  27. %    All Rights Reserved.
  28.     
  29. [nrow,nc]=size(z);nu=th(1,3);ny=nc-nu;
  30. if nargin<5,za=[];end
  31. if nargin<4,nr=[];end
  32. if nargin<3, m=[];end
  33. if isempty(za),za='n';end
  34. if isempty(nr),nr=1:nrow;end
  35. if isempty(m),m=inf;end
  36. nu=th(1,3); if nu==0 & m==inf
  37. error(['For a time series, there must be a finite prediction horizon' ...
  38. '\n(The argument M must be less than inf)']),end
  39. if m==inf,yh=idsim(z(:,ny+1:nc),th);else yh=predict(z,th,m);end
  40. [Ncap,ny]=size(yh);
  41. if za(1)=='y' | za(1)=='Y'
  42. yh(nr,:)=yh(nr,:)-ones(length(nr),1)*yh(nr(1),:);
  43. z(nr,1:ny)=z(nr,1:ny)-ones(length(nr),1)*z(nr(1),1:ny);end
  44. fit=[];
  45. for ky=1:ny
  46. plot([yh(nr,ky) z(nr,ky)])
  47. fittemp = norm(yh(nr,ky)-z(nr,ky))/sqrt(length(nr));    
  48. title(['Output # ',int2str(ky),' Fit: ', num2str(fittemp)])
  49. xlabel('Red/solid: Model output,  Green/dashed: Measured output')
  50. fit=[fit,fittemp];
  51. if ky<ny,pause,end
  52. end
  53.  
  54.  
  55.