Home > wafo > cycles > cmat2extralc.m

cmat2extralc

PURPOSE ^

Extrapolate level crossing spectrum

SYNOPSIS ^

[lcEst,Est,R,MSE] = cmat2extralc(param,F,u,method,plotflag)

DESCRIPTION ^

 CMAT2EXTRALC  Extrapolate level crossing spectrum
 
  CALL: lcEst = cmat2extralc(param,F,u)
        [lcEst,Est,R,MSE] = cmat2extralc(param,F,u,method,plotflag)
 
    param  = Parameter vector, [a b n], defines the grid.
    F      = Cycle matrix. (e.g. Rainflow matrix)          [nxn]
    u      = Levels [u_min u_max], extrapolate below u_min and above u_max. 
    method = A string, describing the method of estimation.
             Generalized Pareto distribution (GPD):
              'gpd,ml'   = Maximum likelihood estimator. (default)
              'gpd,mom'  = Moment method.
              'gpd,pwm'  = Probability Weighted Moments.
             Exponential distribution (GPD with k=0)
              'exp,ml'   = Maximum likelihood estimator.
              'exp,mld'  = Maximum likelihood estimator, discrete.
              'exp,ls1'  = Least Squares estimator, linear, 1 parameter.
              'exp,ls4'  = Least Squares estimator, quadratic, 2 parameters. 
              'exp,wls1' = Weighted Least Squares estimator.
              'exp,wls4' = Weighted Least Squares estimator.
             Rayleigh distribution 
              'ray,ml'   = Maximum likelihood estimator.
    plotflag = 1: Diagnostic plots. (default)
               0: Don't plot diagnostic plots.
 
    lcEst    = The estimated LC.     [struct array]
    Est      = Estimated parameters. [struct array]
    R        = Results from internal calculations. [struct array]
 
  Extrapolates the level crossing spectrum for high and for low levels. 
  The tails of the LC is fited to a survival function of a GPD. 
    H(x) = (1-k*x/s)^(1/k)               (GPD)
  The use of GPD is motivated by POT methods in extreme value theory. 
  For k=0 the GPD is the exponential distribution
    H(x) = exp(-x/s),  k=0               (Exp)
  The tails with the survival function of a Rayleigh distribution.
    H(x) = exp(-((x+x0).^2-x0^2)/s^2)    (Ray)
  where x0 is the value where the LC has its maximum.
  The methods 'gpd,...' uses the GPD. We recommend the use of 'gpd,ml'.
  The methods 'exp,...' uses the Exp. We recommend the use of 'exp,mld'.
  Also the methods 'exp,ls1', 'exp,wls1', 'exp,ls4', 'exp,wls4' work good.
  The method 'ray,ml' uses Ray, and should be used if the load is a Gaussian process.
 
  For further help use 'type cmat2extralc.txt'.
 
  Example: 
    [G,Gh] = mktestmat([-1 1 64],[-0.2 0.2], 0.15,1);
    xD = mctpsim({G Gh},2000);
    Frfc = dtp2rfm(xD,64,'CS');
    [lcEst,Est] = cmat2extralc([-1 1 64],Frfc,[-0.4 0.4]);
    lcG = cmat2lc([-1 1 64],G/sum(sum(G)));
    lcF = cmat2lc([-1 1 64],Frfc);
    clf, semilogx(1000*lcG(:,2),lcG(:,1),lcF(:,2),lcF(:,1),lcEst.lc(:,2),lcEst.lc(:,1))
 
  See also  rfmextrapolate, lc2rfmextreme, extralc, wgpdfit

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

0001 function [lcEst,Est,R,MSE] = cmat2extralc(param,F,u,method,plotflag)
0002 %CMAT2EXTRALC  Extrapolate level crossing spectrum
0003 %
0004 % CALL: lcEst = cmat2extralc(param,F,u)
0005 %       [lcEst,Est,R,MSE] = cmat2extralc(param,F,u,method,plotflag)
0006 %
0007 %   param  = Parameter vector, [a b n], defines the grid.
0008 %   F      = Cycle matrix. (e.g. Rainflow matrix)          [nxn]
0009 %   u      = Levels [u_min u_max], extrapolate below u_min and above u_max. 
0010 %   method = A string, describing the method of estimation.
0011 %            Generalized Pareto distribution (GPD):
0012 %             'gpd,ml'   = Maximum likelihood estimator. (default)
0013 %             'gpd,mom'  = Moment method.
0014 %             'gpd,pwm'  = Probability Weighted Moments.
0015 %            Exponential distribution (GPD with k=0)
0016 %             'exp,ml'   = Maximum likelihood estimator.
0017 %             'exp,mld'  = Maximum likelihood estimator, discrete.
0018 %             'exp,ls1'  = Least Squares estimator, linear, 1 parameter.
0019 %             'exp,ls4'  = Least Squares estimator, quadratic, 2 parameters. 
0020 %             'exp,wls1' = Weighted Least Squares estimator.
0021 %             'exp,wls4' = Weighted Least Squares estimator.
0022 %            Rayleigh distribution 
0023 %             'ray,ml'   = Maximum likelihood estimator.
0024 %   plotflag = 1: Diagnostic plots. (default)
0025 %              0: Don't plot diagnostic plots.
0026 %
0027 %   lcEst    = The estimated LC.     [struct array]
0028 %   Est      = Estimated parameters. [struct array]
0029 %   R        = Results from internal calculations. [struct array]
0030 %
0031 % Extrapolates the level crossing spectrum for high and for low levels. 
0032 % The tails of the LC is fited to a survival function of a GPD. 
0033 %   H(x) = (1-k*x/s)^(1/k)               (GPD)
0034 % The use of GPD is motivated by POT methods in extreme value theory. 
0035 % For k=0 the GPD is the exponential distribution
0036 %   H(x) = exp(-x/s),  k=0               (Exp)
0037 % The tails with the survival function of a Rayleigh distribution.
0038 %   H(x) = exp(-((x+x0).^2-x0^2)/s^2)    (Ray)
0039 % where x0 is the value where the LC has its maximum.
0040 % The methods 'gpd,...' uses the GPD. We recommend the use of 'gpd,ml'.
0041 % The methods 'exp,...' uses the Exp. We recommend the use of 'exp,mld'.
0042 % Also the methods 'exp,ls1', 'exp,wls1', 'exp,ls4', 'exp,wls4' work good.
0043 % The method 'ray,ml' uses Ray, and should be used if the load is a Gaussian process.
0044 %
0045 % For further help use 'type cmat2extralc.txt'.
0046 %
0047 % Example: 
0048 %   [G,Gh] = mktestmat([-1 1 64],[-0.2 0.2], 0.15,1);
0049 %   xD = mctpsim({G Gh},2000);
0050 %   Frfc = dtp2rfm(xD,64,'CS');
0051 %   [lcEst,Est] = cmat2extralc([-1 1 64],Frfc,[-0.4 0.4]);
0052 %   lcG = cmat2lc([-1 1 64],G/sum(sum(G)));
0053 %   lcF = cmat2lc([-1 1 64],Frfc);
0054 %   clf, semilogx(1000*lcG(:,2),lcG(:,1),lcF(:,2),lcF(:,1),lcEst.lc(:,2),lcEst.lc(:,1))
0055 %
0056 % See also  rfmextrapolate, lc2rfmextreme, extralc, wgpdfit
0057 
0058 % References:
0059 %
0060 %   Johannesson, P., and Thomas, J-.J. (2000): 
0061 %   Extrapolation of Rainflow Matrices. 
0062 %   Preprint 2000:82, Mathematical statistics, Chalmers, pp. 18. 
0063 
0064 % Tested  on Matlab  5.3
0065 %
0066 % History:
0067 % Created by PJ (Pär Johannesson) 10-Mar-2000
0068 % Changed by PJ 14-Mar-2000
0069 %   Added sub-function 'make_increasing'
0070 % Changed by PJ 30-Mar-2000
0071 %   Added method 'exp' (GPD with k=0)
0072 % Changed by PJ 05-Apr-2000
0073 %   Modification of 'extrapolation'
0074 % Changed by PJ 11-Apr-2000
0075 %   Added method 'exp,ls'
0076 % Changed by PJ 17-Apr-2000
0077 %   Added methods 'exp,ls1','exp,ls2','exp,ls3','exp,wls1'
0078 % Changed by PJ Apr-2000
0079 %   Added methods 'exp,ls4','exp,wls2','exp,wls3','exp,wls4'
0080 % Changed by PJ 19-Jul-2000
0081 %   Diagnostic plots.
0082 %   Tidy the code.
0083 %   Updated helptext. 
0084 % Updated by PJ 10-Oct-2000
0085 %   More efficient implemantation of 'gpd,ml'.
0086 %   Old version is now called 'gpd,ml0'
0087 %   Some additional changes also made.
0088 % Updated by PJ 03-Nov-2000
0089 %   New method 'ray,ml'
0090 % Updated by PJ 15-Nov-2000
0091 %   Updated method 'gpd,ml'.
0092 %   Now more robust. Does not get into infinite loops.
0093 %   If ML-estimator does not exist, then it uses MOM instead.
0094 % Updated by PJ 6-Dec-2000
0095 %   Updated help text.
0096 % Correction by PJ 11-Dec-2000
0097 %   Now method 'ml' works with Matlab <5.3. ('fzero' changed format)
0098 
0099 % Check input
0100 ni = nargin;
0101 no = nargout;
0102 error(nargchk(3,5,ni));
0103 
0104 if ni < 4, method = []; end
0105 if ni < 5, plotflag = []; end
0106 
0107 % Default values
0108 if isempty(method), method = 'gpd,ml'; end
0109 if isempty(plotflag), plotflag = 1; end
0110 
0111 % Extract level crossings
0112 n = param(3);
0113 paramD = [1 n n];
0114 lc0 = cmat2lc(paramD,F);
0115 uu = levels(param)';
0116 uuD = levels(paramD)';
0117 
0118 % Find the index for the high level
0119 u_high = u(2);
0120 i_high = min(uuD(uu>u_high));
0121 
0122 % Find the index for the low level
0123 u_low = u(1);
0124 i_low = max(uuD(uu<u_low));
0125 
0126 % Remove cycles with  min>i_high  and  max<i_low
0127 F1 = F;
0128 for i = i_high:n
0129   F1(i,:) = 0;
0130 end
0131 for j = 1:i_low
0132   F1(:,j) = 0;
0133 end
0134 
0135 % Get level crossings
0136 lc1 = cmat2lc(paramD,F1);
0137 
0138 lc=lc1;
0139 %lc=lc0;
0140 [dummy,I]=max(lc(:,2));
0141 i_LCmax = lc(I(1),1);
0142 
0143 % Extrapolate LC for high levels
0144 if plotflag, subplot(2,1,1), end
0145 
0146 lcHigh = lc(i_high:end,:);
0147 [lcEst.D.High,Est.D.High,MSE.High] = extrapolate(lcHigh,method,plotflag,i_high-i_LCmax);
0148 
0149 if plotflag, title([method ': Extrapolation for high levels']),  end
0150 
0151 % Extrapolate LC for low levels
0152 if plotflag, subplot(2,1,2), end
0153 
0154 lcLow = lc(1:i_low,:);
0155 lcLow1 = [n+1-flipud(lcLow(:,1)) flipud(lcLow(:,2))];
0156 
0157 [lcEst1,Est1,MSE.Low] = extrapolate(lcLow1,method,plotflag,i_LCmax-i_low);
0158 lcEst.D.Low = [n+1-flipud(lcEst1(:,1)) flipud(lcEst1(:,2))];
0159 Est.D.Low = Est1;
0160 
0161 if isfield(Est.D.Low,'UpperLimit');
0162   Est.D.Low.LowerLimit = n+1-Est.D.Low.UpperLimit;
0163   rmfield(Est.D.Low,'UpperLimit');
0164 end
0165 
0166 if plotflag, title([method ': Extrapolation for low levels']), end
0167 
0168 % Total Level crossing spectrum
0169 lcEst.D.lc = lc;
0170 lcEst.D.lc(1:i_low,:) = lcEst.D.Low;
0171 lcEst.D.lc(i_high:end,:) = lcEst.D.High;
0172 
0173 % Convert to correct scales
0174 lcEst.High = [uu(lcEst.D.High(:,1)) lcEst.D.High(:,2)];
0175 lcEst.Low = [uu(lcEst.D.Low(:,1)) lcEst.D.Low(:,2)];
0176 lcEst.lc = [uu(lcEst.D.lc(:,1)) lcEst.D.lc(:,2)];
0177 
0178 if no > 1
0179   Est.method = method;
0180 end
0181 
0182 % Supplementary results from internal calculations
0183 if no > 2
0184   R.lc0 = cmat2lc(param,F);
0185   R.lc1 = cmat2lc(param,F1);
0186   R.F1 = F1;
0187 end
0188 
0189 % END cmat2extralc
0190 %%%%%
0191 
0192 function [lcEst,Est,MSE] = extrapolate(lc,method,plotflag,offset)
0193 %
0194 % Extrapolate the level crossing specrta for high levels
0195 %
0196   
0197   iu = lc(1,1);
0198   
0199   % Excedences over level u
0200   lc1 = [lc(:,1)-iu lc(:,2)];
0201   lc2 = flipud(lc1);
0202   lc3 = lc2;
0203   
0204   method = lower(method);
0205   methodShape = method(1:3);
0206   methodEst = method(5:end);
0207   
0208   if strcmp(methodShape,'gpd')  & ~strcmp(method,'gpd,ml')
0209   %if strcmp(methodShape,'gpd') | strcmp(method,'exp,ml')| strcmp(method,'exp,mld')
0210     x = [];
0211     for k = 2:length(lc3)
0212       nn = lc3(k,2)-lc3(k-1,2);
0213       if nn > 0
0214         x = [x; lc3(k,1)*ones(nn,1)];
0215       end
0216     end
0217     x = x+.5;   
0218   end
0219   
0220   if strcmp(methodShape,'exp') | strcmp(method,'gpd,ml') | strcmp(method,'ray,ml')
0221     dN = [];
0222     for k = 2:length(lc3)
0223       nn = lc3(k,2)-lc3(k-1,2);
0224       if nn ~= 0
0225         dN = [dN; lc3(k,1) nn]; 
0226       end
0227     end
0228     NN = [dN(:,1) cumsum(dN(:,2))];
0229     NN = flipud(NN);
0230     dN = flipud(dN);
0231   end
0232   
0233   % Estimate tail
0234   if strcmp(methodShape,'gpd') % GPD
0235     
0236     if strcmp(methodEst,'ml')
0237       
0238       dNN = [dN(:,1)+0.5 dN(:,2)];
0239       n = NN(1,2);
0240       
0241       % Calculate start values
0242       X1=dNN(1,1); Xn=dNN(end,1); Xmean = sum(dNN(:,1).*dNN(:,2))/sum(dNN(:,2));
0243       Eps = 1e-6/Xmean;
0244       t_L = 2*(X1-Xmean)/X1^2;     % Lower limit
0245       t_U = 1/Xn-Eps;              % Upper limit
0246       x_L = log(1/Xn - t_L); % Lower limit
0247       x_U = log(1/Xn - t_U); % Upper limit
0248       
0249       x=linspace(x_L,x_U,10);
0250       for i=1:length(x)
0251         f(i)=wgpdfit_mld(x(i),dNN);
0252       end
0253       
0254       I = find(f(1:end-1).*f(2:end) < 0); % Find sign-shift
0255       
0256       % Check if any sign-shift was found
0257       if ~isempty(I) % If ML-estimator exists, then use it! 
0258         i = I(1);
0259         x_start = [x(i) x(i+1)];
0260         if exist('optimset') >= 2 % Function 'optimset' exists ?
0261           % For Matlab 5.3 and higher ???
0262           x_ML = fzero('wgpdfit_mld',x_start,optimset('disp','off'),dNN);
0263         else 
0264           % For Matlab 5.2 and lower ???
0265           x_ML = fzero('wgpdfit_mld',x_start,[],[],dNN);
0266         end
0267         [f,shape,scale] = wgpdfit_mld(x_ML,dNN); % Estimates k_ML and s_ML
0268         
0269         % Calculate the (asymptotic) covariance matrix (obtained by inverting the observed information). 
0270         if (shape >= 1.0),
0271           cov=[];
0272           warning([' The estimate of shape (' num2str(shape) ...
0273               ') is not within the range where the ML estimator is valid (shape<1).'])
0274         elseif (shape >= 0.5),
0275           cov=[];
0276           warning([' The estimate of shape (' num2str(shape) ...
0277               ') is not within the range where the ML estimator is asymptotically normal (shape<1/2).'])
0278         else % Calculate the covariance matrix
0279           Vshape = 1-shape;
0280           Vscale = 2*scale^2;
0281           Covari = scale;
0282           cov = (1-shape)/n*[Vshape Covari; Covari Vscale];
0283         end
0284         
0285       else % If ML-estimator does not exists, then use MOM!
0286         warning([' ML-estimator does not exist for this data. Using MOM instead.'])
0287 
0288         Xvar = sum((dNN(:,1)-Xmean).^2.*dNN(:,2))/(sum(dNN(:,2))-1);
0289         shape = (Xmean^2/Xvar-1)/2;
0290         scale = Xmean*(Xmean^2/Xvar+1)/2;
0291         cov = [];
0292       end
0293       parms = [shape scale];
0294       
0295     else
0296       if strcmp(methodEst,'ml0')
0297         [parms,cov] = wgpdfit(x,'ml',0);
0298       else
0299         [parms,cov] = wgpdfit(x,methodEst,0);
0300       end
0301     end
0302     
0303     Est.k = parms(1);
0304     Est.s = parms(2);
0305     if Est.k>0,
0306       Est.UpperLimit = iu+Est.s/Est.k;
0307     end
0308     Est.cov = cov;
0309     xF = (0:length(lc1)-1)';
0310     F = wgpdcdf(xF,Est.k,Est.s);
0311     
0312   elseif strcmp(methodShape,'exp') % GPD with k=0
0313     
0314     switch methodEst
0315       
0316     case 'ml' % Maximum Likelihood
0317       
0318       % Old version
0319       %n = length(x);
0320       %parms(1) = 0;        % k
0321       %parms(2) = sum(x)/n; % s
0322       
0323       % New version
0324       n = NN(1,2);
0325       Sx = sum(dN(:,2).*(dN(:,1)+0.5));
0326       parms(1) = 0;        % k
0327       parms(2) = Sx/n;     % s
0328       
0329       cov = zeros(2);          % Covariance matrix for estimates.
0330       cov(2,2) = 1/n*parms(2); % Variance of estimate  s
0331       
0332     case 'mld' % Maximum Likelihood, discrete obs.
0333       
0334       n = NN(1,2);
0335       Sx = sum(dN(:,2).*(dN(:,1)+0.0));
0336       s1 = 1/log(1+n/Sx);
0337       
0338       parms(1) = 0;    % k
0339       parms(2) = s1;   % s
0340     
0341       cov = zeros(2);          % Covariance matrix for estimates.
0342       cov(2,2) = 1/n*parms(2); % Variance of estimate  s
0343       
0344       
0345     case {'ls1','wls1'} % One parameter: Least Squares / Weighted LS
0346       
0347       X = NN(:,1);
0348       Y = log(NN(:,2))-log(NN(1,2));
0349       n = length(Y);
0350       
0351       if strcmp(methodEst,'ls1') % Least Squares
0352         W = eye(n);
0353       else % Weighted Least Squares
0354         W = diag(dN(:,2));
0355       end
0356       
0357       % Compute Estimate
0358       th = (X'*W*X)\(X'*W*Y);
0359 
0360       parms(1) = 0;          % k
0361       parms(2) = -1/th(1);   % s
0362       cov = [];
0363             
0364     case {'ls2','wls2'} % Two parameters: Least Squares / Weighted LS
0365           
0366       X = [ones(size(NN(:,1))) NN(:,1)];
0367       Y = log(NN(:,2));
0368       n = length(Y);
0369       
0370       if strcmp(methodEst,'ls2') % Least Squares
0371         W = eye(n);
0372       else % Weighted Least Squares
0373         W = diag(dN(:,2));
0374       end
0375       
0376       % Compute Estimate
0377       th = (X'*W*X)\(X'*W*Y);
0378 
0379       parms(1) = 0;          % k
0380       parms(2) = -1/th(2);   % s
0381       parms(3) = exp(th(1)); % mu0
0382       cov = [];
0383                   
0384     case {'ls3','wls3'} % Three parameters: Least Squares / Weighted LS
0385           
0386       X = [ones(size(NN(:,1))) NN(:,1) NN(:,1).^2];
0387       Y = log(NN(:,2));
0388       n = length(Y);
0389       
0390       if strcmp(methodEst,'ls3') % Least Squares
0391         W = eye(n);
0392       else % Weighted Least Squares
0393         W = diag(dN(:,2));
0394       end
0395       
0396       % Compute Estimate
0397       th = (X'*W*X)\(X'*W*Y);
0398       
0399       parms(1) = 0;          % k
0400       parms(2) = -1/th(2);   % s
0401       parms(3) = exp(th(1)); % mu0
0402       parms(4) = -1/th(3);   % s2
0403       Est.s2 = parms(4);
0404       cov = [];
0405                   
0406     case {'ls4','wls4'} % Three parameters: Least Squares / Weighted LS
0407           
0408       X = [NN(:,1) NN(:,1).^2];
0409       Y = log(NN(:,2))-log(NN(1,2));
0410       n = length(Y);
0411       
0412       if strcmp(methodEst,'ls4') % Least Squares
0413         W = eye(n);
0414       else % Weighted Least Squares
0415         W = diag(dN(:,2));
0416       end
0417       
0418       % Compute Estimate
0419       th = (X'*W*X)\(X'*W*Y);
0420       
0421       parms(1) = 0;          % k
0422       parms(2) = -1/th(1);   % s
0423       parms(4) = -1/th(2);   % s2
0424       Est.s2 = parms(4);
0425       cov = [];
0426             
0427     end % switch methodEst
0428     
0429     Est.k = parms(1);
0430     Est.s = parms(2);
0431     Est.cov = cov;
0432     xF = (0:length(lc1)-1)';
0433     switch methodEst
0434     case {'ml', 'mld', 'ls1', 'wls1', 'ls2', 'wls2'}
0435       F = 1 - exp(-xF/Est.s);
0436     case {'ls3', 'wls3', 'ls4', 'wls4'}
0437       F = 1 - exp(-xF/Est.s-xF.^2/Est.s2);
0438     end
0439     
0440   elseif strcmp(methodShape,'ray') % Rayleigh
0441     
0442     n = NN(1,2);
0443     Sx = sum(dN(:,2).*((dN(:,1)+offset+0.5).^2-(offset)^2));
0444     a=sqrt(Sx/n);  % Shape parameter
0445     
0446     Est.s = a;
0447     
0448     xF = (0:length(lc1)-1)';
0449     F = 1 - exp(-((xF+offset).^2-offset^2)/a^2);
0450     
0451   else % Unknown method
0452     
0453     error(['Unknown method: ' method]);
0454     
0455   end
0456   
0457   % Estimate of mu0
0458   switch methodEst
0459   case {'ls2','wls2','ls3','wls3'}
0460     Est.lcu = parms(3);
0461   otherwise
0462     Est.lcu = lc(1,2);
0463   end
0464   
0465   lcEst = [xF+iu Est.lcu*(1-F)];
0466   lcEst(end,2) = 0;  % No crossings of the highest level
0467   
0468   % Compute MSE (Mean Square Error)
0469   if strcmp(method,'gpd,ml') | strcmp(method,'exp,ml') | strcmp(method,'exp,mld') | strcmp(method,'ray,mld')
0470     MSE = 0;
0471     % Not impemented
0472   elseif strcmp(methodShape,'gpd')
0473     n=length(x);
0474     q1 = ((1:n)-1/2)/n;
0475     xq1 = wgpdinv(q1,Est.k,Est.s)';
0476     MSE = sum((x-xq1).^2)/n;
0477   elseif strcmp(method,'ray,ml')
0478     MSE = 0;
0479     % Not impemented
0480   else % (Weighted) Least Squares estimate
0481     Lam = sqrt(W);
0482     MSE = (Lam*X*th - Lam*Y)'*(Lam*X*th - Lam*Y)/sum(diag(W));
0483   end
0484   
0485   % Plot Diagnostics 
0486   if plotflag
0487     if strcmp(method,'gpd,ml') | strcmp(method,'exp,ml')| strcmp(method,'exp,mld')
0488       Csum = cumsum(dN(:,2));
0489       q1 = ([1; Csum(1:end-1)+1; Csum(1:end)]-1/2)/n;
0490       xq1 = wgpdinv(q1,Est.k,Est.s)';
0491       x = [dN(:,1); dN(:,1)]+0.5;
0492       wqqplot(x,xq1); hold on
0493       plot([0 max(x)],[0 max(x)]), hold off
0494       xlabel('data'), ylabel('Quantile')
0495     elseif strcmp(methodShape,'gpd')
0496       wqqplot(x,xq1); hold on
0497       plot([0 max(x)],[0 max(x)]), hold off
0498       xlabel('data'), ylabel('Quantile')
0499     elseif strcmp(method,'ray,ml')
0500       Csum = cumsum(dN(:,2));
0501       q1 = ([1; Csum(1:end-1)+1; Csum(1:end)]-1/2)/n;
0502       xq1 = sqrt(offset^2-Est.s^2*log(1-q1))-offset;
0503       x = [dN(:,1); dN(:,1)]+0.5;
0504       wqqplot(x,xq1); hold on
0505       plot([0 max(x)],[0 max(x)]), hold off
0506       xlabel('data'), ylabel('Quantile')
0507     else
0508       plot(X(:,1),Y,'+'), hold on
0509       plot(X(:,1),X*th,'r'), hold off
0510       xlabel('data'), ylabel('log(lc)')
0511     end
0512   end
0513     
0514 % END extrapolate
0515 %%%%%
0516

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