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:
 wgpdcdf Generalized Pareto cumulative distribution function wgpdfit Parameter estimates for Generalized Pareto data eig Eigenvalues and eigenvectors. error Display message and abort function. interp1 1-D interpolation (table lookup) linspace Linearly spaced vector. lower Convert string to lowercase. mean Average or mean value. nan Not-a-Number. semilogx Semi-log scale plot.
This function is called by:
 test_cycles Quick test of the routines in module 'cycles'

## 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