Home > wafo > cycles > extralc.m

extralc

PURPOSE ^

Extrapolate level crossing spectrum

SYNOPSIS ^

[lcEst,Est] = extralc(lc,u,method,plotflag)

DESCRIPTION ^

 EXTRALC  Extrapolate level crossing spectrum
 
  CALL: [lcEst,Est] = extralc(lc,u,method,plotflag)
 
    lc     = Level crossing spectrum.                       [n,2]
    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.
              'gpd,pkd'  = Pickands' estimator.
             Exponential distribution (GPD with k=0)
              'exp,ml'   = Maximum likelihood 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]
 
  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 method 'exp,ml' uses the Exp. 
  The method 'ray,ml' uses Ray, and should be used if the load is a Gaussian process.
 
  Example: 
    S = jonswap;
    x = spec2sdat(S,100000,0.1,[],'random');
    lc = dat2lc(x); s = std(x(:,2));
    [lcEst,Est] = extralc(lc,s*[-2 2]);
    [lcEst,Est] = extralc(lc,s*[-2 2],'exp,ml');
    [lcEst,Est] = extralc(lc,s*[-2 2],'ray,ml');
 
  See also  cmat2extralc, rfmextrapolate, lc2rfmextreme, extralc, wgpdfit

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

001 function [lcEst,Est] = extralc(lc,u,method,plotflag)
002 %EXTRALC  Extrapolate level crossing spectrum
003 %
004 % CALL: [lcEst,Est] = extralc(lc,u,method,plotflag)
005 %
006 %   lc     = Level crossing spectrum.                       [n,2]
007 %   u      = Levels [u_min u_max], extrapolate below u_min and above u_max. 
008 %   method = A string, describing the method of estimation.
009 %            Generalized Pareto distribution (GPD):
010 %             'gpd,ml'   = Maximum likelihood estimator. (default)
011 %             'gpd,mom'  = Moment method.
012 %             'gpd,pwm'  = Probability Weighted Moments.
013 %             'gpd,pkd'  = Pickands' estimator.
014 %            Exponential distribution (GPD with k=0)
015 %             'exp,ml'   = Maximum likelihood estimator.
016 %            Rayleigh distribution 
017 %             'ray,ml'   = Maximum likelihood estimator.
018 %   plotflag = 1: Diagnostic plots. (default)
019 %              0: Don't plot diagnostic plots.
020 %
021 %   lcEst    = The estimated LC.     [struct array]
022 %   Est      = Estimated parameters. [struct array]
023 %
024 % Extrapolates the level crossing spectrum for high and for low levels. 
025 % The tails of the LC is fited to a survival function of a GPD. 
026 %   H(x) = (1-k*x/s)^(1/k)               (GPD)
027 % The use of GPD is motivated by POT methods in extreme value theory. 
028 % For k=0 the GPD is the exponential distribution
029 %   H(x) = exp(-x/s),  k=0               (Exp)
030 % The tails with the survival function of a Rayleigh distribution.
031 %   H(x) = exp(-((x+x0).^2-x0^2)/s^2)    (Ray)
032 % where x0 is the value where the LC has its maximum.
033 % The methods 'gpd,...' uses the GPD. We recommend the use of 'gpd,ml'.
034 % The method 'exp,ml' uses the Exp. 
035 % The method 'ray,ml' uses Ray, and should be used if the load is a Gaussian process.
036 %
037 % Example: 
038 %   S = jonswap;
039 %   x = spec2sdat(S,100000,0.1,[],'random');
040 %   lc = dat2lc(x); s = std(x(:,2));
041 %   [lcEst,Est] = extralc(lc,s*[-2 2]);
042 %   [lcEst,Est] = extralc(lc,s*[-2 2],'exp,ml');
043 %   [lcEst,Est] = extralc(lc,s*[-2 2],'ray,ml');
044 %
045 % See also  cmat2extralc, rfmextrapolate, lc2rfmextreme, extralc, wgpdfit
046 
047 % References:
048 %
049 %   Johannesson, P., and Thomas, J-.J. (2000): 
050 %   Extrapolation of Rainflow Matrices. 
051 %   Preprint 2000:82, Mathematical statistics, Chalmers, pp. 18. 
052 
053 % Tested  on Matlab  5.3
054 %
055 % History:
056 % Created by PJ (Pär Johannesson) 10-Mar-2000
057 % Changed by PJ 14-Mar-2000
058 %   Added sub-function 'make_increasing'
059 % Changed by PJ 15-Mar-2000
060 %   Added sub-function 'make_increasing'
061 % Updated by PJ 07-Dec-2000
062 %   Added 'exp,ml' and 'ray,ml'
063 %   Help text
064 % Corrected by PJ 17-Feb-2004
065 
066 % Check input
067 
068 ni = nargin;
069 no = nargout;
070 error(nargchk(1,4,ni));
071 
072 if ni < 3, method = []; end
073 if ni < 4, plotflag = []; end
074 
075 % Default values
076 if isempty(method), method = 'gpd,ml'; end
077 if isempty(plotflag), plotflag = 1; end
078 
079 % Maximum of lc
080 [M,I] = max(lc(:,2));  % Corrected by PJ 17-Feb-2004
081 lc_max = lc(I(1),1);
082 
083 % Extrapolate LC for high levels
084 [lcEst.High,Est.High] = extrapolate(lc,u(2),method,u(2)-lc_max);
085 
086 % Extrapolate LC for low levels
087 
088 lc1 = [-flipud(lc(:,1)) flipud(lc(:,2))];
089 
090 [lcEst1,Est1] = extrapolate(lc1,-u(1),method,lc_max-u(1));
091 lcEst.Low = [-flipud(lcEst1(:,1)) flipud(lcEst1(:,2:end))];
092 Est.Low = Est1;
093 
094 if plotflag
095   semilogx(lc(:,2),lc(:,1),lcEst.High(:,2),lcEst.High(:,1),lcEst.Low(:,2),lcEst.Low(:,1))
096 end
097 
098 %%%
099 function [lcEst,Est] = extrapolate(lc,u,method,offset)
100   % Extrapolate the level crossing specrta for high levels
101   
102   method = lower(method);
103   methodShape = method(1:3);
104   methodEst = method(5:end);
105   
106   % Excedences over level u
107   Iu = lc(:,1)>u;
108   lc1 = lc(Iu,:);
109   lc2 = flipud(lc1);
110   lc3 = make_increasing(lc2);
111   
112   % Corrected by PJ 17-Feb-2004
113   lc3=[0 0; lc3]; x=[];
114   for k=2:length(lc3(:,1))
115       nk = lc3(k,2)-lc3(k-1,2);
116       x = [x; ones(nk,1)*lc3(k,1)];
117   end
118   x = x-u;   
119   
120   % Estimate tail
121   switch methodShape
122     
123   case 'gpd'
124     
125     [parms,cov] = wgpdfit(x,methodEst,0);
126     
127     Est.k = parms(1);
128     Est.s = parms(2);
129     Est.cov = cov;
130     if Est.k>0,
131       Est.UpperLimit = u+Est.s/Est.k;
132     end
133     
134     xF = (0.0:0.01:4)';
135     F = wgpdcdf(xF,Est.k,Est.s);
136     Est.lcu = interp1(lc(:,1),lc(:,2),u)+1;
137     
138     % Calculate 90 % confidence region, an ellipse, for (k,s)
139     [B,D] =eig(cov);
140     b = [Est.k; Est.s];
141     
142     r = sqrt(-2*log(1-90/100)); % 90 % confidence sphere
143     Nc = 16+1;
144     ang = linspace(0,2*pi,Nc);
145     c0 = [r*sqrt(D(1,1))*sin(ang); r*sqrt(D(2,2))*cos(ang)]; % 90% Circle
146 %    plot(c0(1,:),c0(2,:))
147     
148     c1 = B*c0+b*ones(1,length(c0)); % Transform to ellipse for (k,s)
149 %    plot(c1(1,:),c1(2,:)), hold on
150     
151     % Calculate conf.int for lcu
152     % Assumtion: lcu is Poisson distributed
153     % Poissin distr. approximated by normal when calculating conf. int.
154     dXX = 1.64*sqrt(Est.lcu); % 90 % quantile for lcu
155     
156     lcEstCu = zeros(length(xF),Nc); lcEstCl = lcEstCu;
157     for i=1:Nc
158       k=c1(1,i);
159       s=c1(2,i);
160       F2 = wgpdcdf(xF,k,s);
161       lcEstCu(:,i) = (Est.lcu+dXX)*(1-F2);
162       lcEstCl(:,i) = (Est.lcu-dXX)*(1-F2);
163     end
164     
165     lcEstCI = [min(lcEstCl')' max(lcEstCu')'];
166     
167     lcEst = [xF+u Est.lcu*(1-F) lcEstCI];
168     
169   case 'exp'
170     
171     n = length(x);
172     s = mean(x);
173     cov = s/n;
174     Est.s = s;
175     Est.cov = cov;
176     
177     xF = (0.0:0.01:4)';
178     F = 1-exp(-xF/s);
179     Est.lcu = interp1(lc(:,1),lc(:,2),u)+1;
180     
181     lcEst = [xF+u Est.lcu*(1-F)];
182     
183   case 'ray'
184     
185     n = length(x);
186     Sx = sum((x+offset).^2-offset^2);
187     s=sqrt(Sx/n);  % Shape parameter
188     
189     Est.s = s;
190     Est.cov = NaN;
191     
192     xF = (0.0:0.01:4)';
193     F = 1 - exp(-((xF+offset).^2-offset^2)/s^2);
194     Est.lcu = interp1(lc(:,1),lc(:,2),u)+1;
195     
196     lcEst = [xF+u Est.lcu*(1-F)];
197     
198   otherwise
199     
200     error(['Unknown method: ' method]);
201     
202   end
203   
204     
205 %% End extrapolate
206   
207 function y=make_increasing(x);
208   % Makes the signal x strictly increasing. 
209   %
210   % x = two column matrix with times in first column and values in the second.
211   
212   n = length(x);
213   
214   i=1;
215   y = x;
216   y(1,:) = x(1,:); j=1;
217   
218   while i<n
219     while x(i,2)<=y(j,2) & i<n
220       i = i+1;
221     end
222     if x(i,2)>y(j,2)
223       j = j+1;
224       y(j,:) = x(i,:);
225     end
226   end
227   
228   y = y(1:j,:);
229   
230 %% End make_increasing
231 
232

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