Home > wafo > damage > snplot.m

snplot

PURPOSE ^

Plots SN-data and estimates parameters

SYNOPSIS ^

[e,m,sigma2]= snplot(S,N,nr)

DESCRIPTION ^

 SNPLOT Plots SN-data and estimates parameters 
  according to the least-squares method
 
  CALL:  [e,beta,s2] = snplot(S,N,nr);
 
   where
 
         e,beta = the parameters in the model, 
         s2     = the residual variance for linear regression 
                  in logS/logN plane,
         S      = a nx1 vector with S-data,
         N      = a nx1 vector with N-data,
         nr     = plot parameter (optional input argument);
 
                  1 = only SN-data will be plotted,
                  2 = SN-data and fitted line will be plotted,
                  3 = only log(S)/log(N)-data will be plotted,
                  4 = log(S)/log(N)-data and fitted line will be plotted,
                 11-14 = Same as above but x-axis = N, y-axis = S . 
   
  Model:
        N(s) = K/(e*s^beta)
 
  Example:  
    sn = load('sn.dat'); s = sn(:,1); N = sn(:,2);
    [e,beta,s2] = snplot(s,N,4)   % S-N, x-axis = S, y-axis = N
    [e,beta,s2] = snplot(s,N,14)  % N-S, x-axis = N, y-axis = S

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function [e,m,sigma2]= snplot(S,N,nr)
002 %SNPLOT Plots SN-data and estimates parameters 
003 % according to the least-squares method
004 %
005 % CALL:  [e,beta,s2] = snplot(S,N,nr);
006 %
007 %  where
008 %
009 %        e,beta = the parameters in the model, 
010 %        s2     = the residual variance for linear regression 
011 %                 in logS/logN plane,
012 %        S      = a nx1 vector with S-data,
013 %        N      = a nx1 vector with N-data,
014 %        nr     = plot parameter (optional input argument);
015 %
016 %                 1 = only SN-data will be plotted,
017 %                 2 = SN-data and fitted line will be plotted,
018 %                 3 = only log(S)/log(N)-data will be plotted,
019 %                 4 = log(S)/log(N)-data and fitted line will be plotted,
020 %                11-14 = Same as above but x-axis = N, y-axis = S . 
021 %  
022 % Model:
023 %       N(s) = K/(e*s^beta)
024 %
025 % Example:  
026 %   sn = load('sn.dat'); s = sn(:,1); N = sn(:,2);
027 %   [e,beta,s2] = snplot(s,N,4)   % S-N, x-axis = S, y-axis = N
028 %   [e,beta,s2] = snplot(s,N,14)  % N-S, x-axis = N, y-axis = S
029 
030   
031 % Tested on: Matlab 6.0
032 % History:
033 % Correction by PJ 07-Jul-2005
034 %   Changed 'break' to 'return'
035 % Revised by jr 01-Apr-2001
036 % - example, help text
037 % Revised by PJ 10-Jan-2000
038 %   updated for WAFO
039 % Original version from FAT by Mats Frendahl 
040 %   Copyright 1993, Mats Frendahl, Dept. of Math. Stat., University of Lund.
041 
042 [n,m]=size(S);
043 if m>n, S=S'; end
044 [n,m]=size(S);
045 if m~=1, disp('   First argument must be  nx1, program will terminate.'), end
046 
047 [n,m]=size(N);
048 if m>n, N=N'; end
049 [n,m]=size(N);
050 if m~=1, disp('   Second argument must be  nx1, program will terminate.'), end
051 
052 A=[ones(length(S),1) log(S)];
053 beta=inv(A'*A)*A'*log(N); m=-beta(2); e=exp(-beta(1));
054 Q0=log(N)'*log(N)-beta'*A'*log(N);
055 sigma2=Q0/(length(log(S))-2);
056 k=exp(sigma2/2);
057 borderS=[min(S) max(S)]; 
058 borderN=[min(N) max(N)]; 
059 
060 number_of_s=99;
061 
062 if nargin==3
063   if nr<10
064     if nr==1 | nr==3
065      plot(S,N,'.','markersize',12)
066      title('SN-data')
067     elseif nr==2 | nr==4
068       s=borderS(1):(borderS(2)-borderS(1))/number_of_s:borderS(2);
069       n=k/e*s.^(-m);
070       plot(s,n,'-',S,N,'.','markersize',12)
071       axis([[0.9 1.1].*borderS 0 1.1*borderN(2)]);
072       title('SN-data with estimated N(s)'),xlabel('s'), ylabel('N(s)')
073     end
074     
075     axis([[0.9 1.1].*borderS 0 1.1*borderN(2)]);
076     xlabel('s'), ylabel('N')
077     
078   else
079     if nr==11 | nr==13
080       plot(N,S,'.','markersize',12)
081       title('SN-data')
082     elseif nr==12 | nr==14
083       s=borderS(1):(borderS(2)-borderS(1))/number_of_s:borderS(2);
084       n=k/e*s.^(-m);
085       plot(n,s,'-',N,S,'.','markersize',12)
086       axis([[0.9 1.1].*borderS 0 1.1*borderN(2)]);
087       title('SN-data with estimated N(s)')
088     end
089     
090     axis([0 1.1*borderN(2) [0.9 1.1].*borderS ]);
091     xlabel('N'), ylabel('s')
092   end
093   
094   if mod(nr,10) > 2 % if number == 3,4,13,14, then log-scale
095     h=gca;
096     set(h,'YScale','log','XScale','log');
097   end
098   
099 end
100 
101 return
102 
103 %
104 % Old plots
105 %
106 if nargin==3
107    if nr==1
108      plot(S,N,'.','markersize',12)
109      axis([[0.9 1.1].*borderS 0 1.1*borderN(2)]);
110      title('SN-data'),xlabel('s'), ylabel('N')
111    elseif nr==2
112      s=borderS(1):(borderS(2)-borderS(1))/number_of_s:borderS(2);
113      n=k/e*s.^(-m);
114      plot(s,n,'-',S,N,'.','markersize',12)
115      axis([[0.9 1.1].*borderS 0 1.1*borderN(2)]);
116      title('SN-data with estimated N(s)'),xlabel('s'), ylabel('N(s)')
117    elseif nr==3
118      axis([[0.9 1.1].*log(borderS) 0 1.1*log(borderN(2))]);
119      plot(log(S),log(N),'.','markersize',12)
120      title('log(S)/log(N)-data'),xlabel('log(s)'), ylabel('log(N(s))')
121    elseif nr==4
122      y=beta(1)+beta(2)*log(borderS);
123      axis([[0.9 1.1].*log(borderS) 0 1.1*log(borderN(2))]);
124      plot(log(borderS),y,'-',log(S),log(N),'.','markersize',12)
125      title('log(S)/log(N)-data with estimated log(N(s))')
126      xlabel('log(s)'), ylabel('log(N(s))')
127    end
128 end
129 
130

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003